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PREFACE 


This  report  is  the  fifth  in  a  series  of  reports  on  constitutive  modeling  by 
Applied  Research  Associates,  Inc.  (ARA)  which  have  been  funded  by  the  Air  Force 
Office  of  Scientific  Research  (AFOSR).  The  research  described  by  these  reports 
and  other  publications  has  been  directed  toward  improved  calculational  modeling 
of  soils  under  complex  dynamic  loadings,  specifically  those  produced  by 
explosive  sources.  The  following  sunrnary  of  these  ARA  reports  is  intended  to 
provide  background  and  perspective  for  the  current  report. 

The  first  report  [Oass,  Bratton,  and  Higgins  (1981)]  was  primarily  a  review 
of  constitutive  modeling  requirements  for  dynamic  modeling  of  soil  behavior.  A 
literature  review  of  existing  models  was  presented  to  show  which  models  or  parts 
of  models  might  be  applied  to  the  specific  problem  of  explosive  loadings.  The 
Soil  Element  Model  (SEM),  a  utility  computer  program  used  to  study  .and  develop 
material  models,  was  also  introduced.  The  second  report  [Dass,  Merkle,  and 
Bratton  (1983)]  dealt  with  the  capability  of  several  selected  models  to  predict 
the  behavior  of  soils  under  one-dimensional  planar,  cylindrical,  and  spherical 
geometry  explosive  loadings.  The  important  behavioral  differences  between 
models  whose  parameters  were  fit  to  laboratory  data  and  those  fit  to  insitu  data 
were  illustrated.  The  third  report  [Merkle  and  Dass  (1983)]  focused  on  modeling 
the  dynamic  response  of  saturated  soil  and  a  review  of  some  widely-applied 
plasticity  concepts.  It  provided  substantial  theoretical  background  toward 
development  of  an  Improved  constitutive  model,  described  In  the  fourth  report 
[Merkle  and  Dass  (1985)].  The  new  model  was  based  on  work  by  Lade  at  UCLA  with 
improvements  aimed  at  better  response  in  a  single-phase  finite  difference 
calculation  of  explosi vely-dri ven  wave  propagation.  In  the  fourth  report,  the 
detailed  theoretical  background  of  this  model  and  important  aspects  of  several 
other  models  were  presented,  and  its  single-element  behavior  was  directly 
compared,  using  the  SEM,  to  other  models  currently  in  use  in  the  ground  shock 
corrmunity. 

This  report  describes  the  first  calculational  tests  of  the  new  model  (now 
referred  to  as  the  ARA  Three-Invariant  Model)  for  one-  and  two-dimensional  wave 
propagation.  It  addresses  many  of  the  computational  issues  which  need  to  be 
considered  when  fully  implementing  a  constitutive  model.  The  theory  behind  the 
model  is  not  presented  in  detail  and  the  reader  may  find  it  necessary  to  consult 
the  previous  reports,  particularly  the  fourth  report,  for  additional  background 
information. 
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1.0  INTRODUCTION 


The  area  of  constitutive  modeling  in  soil  mechanics  and  dynamics  has 
evolved  to  the  point  where  many  kinds  of  models  have  been  proposed,  fewer  have 
been  implemented,  and  still  fewer  are  actually  being  used.  There  are  several 
reasons  for  this: 

(1)  Models  which  have  been  around  for  a  while  are  used  more  often  because 
they  are  familiar,  ready  to  use,  their  limitations  are  known,  and 
organizations  have  developed  an  experience  base  over  the  course  of 
many  projects. 

(2)  Complicated  models  are  harder  to  understand  and  therefore  it  is  often 
hard  to  attach  physical  significance  to  their  features  and 
parameters. 

(3)  When  implemented,  complex  models  sometimes  do  not  produce  appreciably 
better  calculated  results. 

(4)  Involved  laboratory  testing  for  determining  model  parameters  is  beyond 
the  scope  of  many  projects. 

By  way  of  analogy,  the  selection  of  a  constitutive  model  for  use  in  an  engineer¬ 
ing  calculation  can  be  compared  to  the  purchase  of  an  automobile.  The  engineer 
is  faced  with  choices,  many  of  which  are  similar  to  those  faced  by  the  car-buyer 
(see  Table  1).  The  actual  choice  usually  boils  down  to  budget,  prior  exper¬ 
ience,  or  the  advice  of  a  helpful  friend.  For  a  model  to  gain  acceptance  in  the 
marketplace,  the  model  developer  must  be  aware  of  what  model-users  are  looking 
for  and  what  they  need. 

The  purpose  of  this  research  has  been  to  test  a  material  model  (developed 
under  a  prior  AFOSR  contract)  in  much  the  same  fashion  as  a  car  manufacturer 
would  test  a  new  model  before  making  it  available  to  the  consumer.  The  model 
has  been  taken  over  some  bumpy  calculational  roads  to  be  sure  the  doors  won’t 
fall  off  in  the  process.  Its  usability  and  speed  have  been  tested  and  compared 
with  some  other  models. 

The  result  is  a  model  with  no  known  bugs,  some  improvements,  and  initial 
computational  experience.  The  model's  future  now  depends  to  a  large  degree  on 
Its  attractiveness  to  potential  consumers  (i.e.,  calculators).  Many  of  its 
features  have  never,  to  the  authors'  knowledge,  actually  been  implemented  in 
calculations  of  blast  effects.  There  is  promise  that  this  model  will  produce 
better  calculations  of  dynamic  problems  involving  complex  loading  paths. 


TABLE  1.  PARALLELS  IN  MODEL  EVALUATION 


Choosing  An 

Automobile1 

Choosing  A 

Constitutive  Model 

Performance 

Fuel  Economy 
Driveability 
Acceleration 

Braking 

Hand! 1 ng/Roadhol dl ng 

Computer  Funds  Econorry 
Useablllty 

Speed 

Theoretical  Soundness 
Predictive  Accuracy 

Performance 

Driver 

Driving  Position 

I nstruments/Control s 

Visibility 

Heating/Ventilation 

Installed  In  Right  Code 
Parameter  Determination 
Physical  Insight  Into  Model 
Source  for  Help  In  Problems 

User 

Passenger 

Seat  Comfort 

Passenger  Room 

Ride 

Noise 

Confidence  In  Model 

Physical  Insight  Into  Model 
Quality  of  Results 

Model  Induced  Errors 

Client 

Convenience 

Entry /Ex it 

Cargo  Room 

Serviceability 

Equipment 

Learning  Time 

Expandability 

Serviceability 

Available  Behavioral 

Features 

Convenience 

Workmanship 

Body  Construction 

Paint/Exterior 

Interior 

Model  Construction 
Implementation 

Coding  Slmpl 1  city /Accuracy 

Workmanship 

1  Consumer  Guide:  1985  Cars,  Publications  International,  Ltd.,  1985,  p.  376. 
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2.0  ARA  THREE  INVARIANT  MODEL 


2.1  Model  Background 


The  ARA  three  invariant  model  was  developed  by  Applied  Research  Associates, 
Inc.  (ARA)  under  funding  from  the  Air  Force  Office  of  Scientific  Research 
(AFOSR).  The  model  Is  described  In  detail  and  compared  with  several  other 
models  under  laboratory  loading  conditions  by  Merkle  and  Dass  (1985).  The  model 
owes  Its  beginnings  and  concept  to  a  plasticity  model  developed  by  Lade  and 
Nelson  (1981)  at  UCLA.  The  ARA  model  is  called  a  three  Invariant  model  because 
the  controlling  surfaces  in  principal  stress  space  are  functions  of  three 
independent  stress  Invariants  (octahedral  normal  and  shear  stress  and  Lode's 
angle).  Figure  1  shows  several  views  of  the  model’s  surfaces. 


The  ARA  model  is  a  strain  hardening/softening  elastoplastic  model  with  two 
independent  yield  surfaces,  one  associative  and  the  other  nonassociati ve.  The 
compressive  yield  surface  Is  an  ellipsoid  with  its  center  at  the  origin  In 
principal  stress  space.  It  is  associative,  only  strain  hardens,  and  the  strain 
hardening  parameter  Is  the  corresponding  plastic  work.  The  expansive  yield 
surface  is  a  hyperboloid  with  its  apex  on  the  hydrostatic  axis  in  principal 
stress  space.  It  Is  nonassociati ve,  both  strain  hardens  and  softens,  and  the 
strain  hardening/softening  parameter  is  the  corresponding  plastic  work.  The 
expansive  plastic  potential  surface  is  also  a  hyperboloid. 


The  ARA  model  computational  strategy  involves  four  Independent  conditions 
for  each  yield  surface: 


(1)  Yield  Condition  -  This  is  a  necessary,  but  not  sufficient  condition 
for  yielding  to  occur. 


(2)  Flow  Rule  -  If  yielding  occurs,  the  plastic  strain  Increment  is  normal 
to  the  plastic  potential  surface. 


(3)  Consistency  Condition  -  If  yielding  occurs,  the  yield  condition  must 
be  satisfied  throughout  yielding. 

(4)  Dissipation  Condition  -  Yielding  must  generate  positive  plastic  work. 


The  basic  equations  of  the  ARA  model  without  tensile  strength  are  given 
below.  The  compressive  yield  criterion  is 


fc  *  fi 


-  f"  =  0 

C 


where 
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Principal  Stress  Space  Perspective. 


Figure  1.  ARA  Three-Invariant  Model  Yield 
and  Plastic  Potential  Surfaces. 


fc  =  pa 


2  /  "c  \ 
VC  P  J 


where  f'  =  compressive  yield  stress  function,  f”  =  compressive  yield  hardening 
function,  aoct  =  octahedral  normal  stress,  ioct  =  octahedral  shear  stress, 

Pa  =  atmospheric  pressure,  Wc  =  compressive  plastic  work  (defined  below),  and 
r,  C,  P  =  model  parameters.  The  compressive  flow  rule  Is 


<dec)  =  d\c  j~| 


where  (d tc}  =  column  vector  of  compressive  plastic  strain  Increments, 
dXc  =  compressive  yield  proportional ity  constant,  (<j)  =  column  vector  of 
effective  stress  components.  The  compressive  plastic  work  Increment  is 

dWc  =  (c)T  (dtc)  >  0 

The  expansive  yield  criterion  is 

*  =  f‘  -  f"  =  0 

P  P  P 


where 


✓Toct\  /  pa  \ 

f-(— )(l-£c°s  3.)(—  ♦  ■) 


u_  1/q 

p  Vp„/ 


f "  =  n 

p,max  nl 


where  w  =  Lode's  angle,  and  E,  m,  n^,  a,  b,  q  =  model  parameters.  The  expansive 
flow  rule  Is 
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where  {dtp}  =  column  vector  of  expansive  plastic  strain  increments, 

d\p  =  expansive  yield  proportionality  constant.  The  expansive  plastic  potential 

function  Is 


'°OCt\ 


9p 


/toct\ 

=  ( — — )  (1  -  E  cos  3u)  - 


V  Pa  > 


(U) 


,  /°oct\ 
1  +  m(— ) 


where  =  model  parameter.  The  expansive  plastic  work  increment  Is 

dWp  =  {dtp}  >  0 


(12) 


The  effective  stress  increments  are  determined  by  the  elastic  strain 
increments 


(do)  =  Ce  (dee) 


(13) 


where 


(dee)  =  (de)  -  (dec)  -  (dcp)  (14) 

where  Ce  =  elastic  stiffness  matrix,  (dee)  =  column  vector  of  elastic  strain 
increments,  and  (de)  =  column  vector  of  total  strain  increments. 


2.2  Low  Stress  Improvements 

During  the  course  of  the  continuum  code  checkout,  the  ARA  model  was 
improved  in  several  ways.  One  of  the  Improvements  is  an  option  for  tensile 
strength  (cohesion).  Originally,  the  expansive  failure  surface  had  its  apex  at 
the  origin  In  principal  stress  space.  This  precluded  modeling  tensile  strength, 
which  Is  exhibited  by  some  soils  and  most  rocks.  Tension  capacity  was  added  by 
shifting  the  entire  model  down  the  hydrostatic  axis,  as  shown  in  Figure  2. 

The  new  equation  for  the  expansive  yield  surface  is  then 

^oct\  _  _  /Pa 

3a 


%  ’  (-f1)  u  - E  cos  3“)  ( 


°oct+T 


+  m 


(15) 


where  T  Is  tensile  strength  measured  along  the  octahedral  normal  stress  axis. 
The  revised  equation  for  the  expansive  plastic  potential  is 


The  compressive  yield  surface  is  also  shifted: 

fc  -  3< vt  +  T>2  +  3  r2  Vt2  (H) 

The  net  effect  Is  a  simple  translation  of  the  principal  stress  axes.  Equations 
which  Involve  the  derivatives  of  the  above  quantities  are  also  affected  and  have 
been  modified. 

Figure  3  compares  the  response  of  the  model  using  zero  tension  capacity 
with  that  using  2.0  MPa  tension  capacity  for  an  element  exercised  along  an 
arbitrary  strain  path.  The  strain  path  Is  representative  of  those  experienced 
in  spherical  wave  propagation  [Akers  (1985)]. 

When  soil  is  subjected  to  large  tensile  strains,  the  soil  particles 
separate  and  the  material  behaves  less  and  less  like  a  continuum  as  large  voids 
develop.  This  kind  of  behavior  is  modeled  in  the  ARA  model  by  tracking  the 
total  volume  strain  which  occurs  while  the  element  is  failed  in  tension.  An 
element  is  not  allowed  to  rejoin  (develop  compressive  stress)  until  the  volume¬ 
tric  strain  is  equal  to  that  at  which  tensile  failure  occurred.  Specifically, 
the  strain  tracked  is 

e  spa 1 1  =  cj  +  c2  +  c 3  (18) 

While  the  material  is  in  a  spalled  condition,  each  principal  stress  is  set  equal 
to  -T,  and  all  shear  stresses  are  set  equal  to  zero.  Thus,  the  stress  point 
remains  at  the  apex  of  the  expansive  yield  surface  and  rejoin  always  initiates 
from  this  point.  Figure  4  shows  the  behavior  of  an  element  which  has  failed  in 
tension  (using  the  same  strain  path  shown  in  Figure  3)  and  is  then  forced  to 
rejoin  under  subsequent  uniaxial  loading  (ej  =  C3  =  0,  -  compressi ve) . 
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Another  substantial  improvement  to  the  ARA  model  is  the  addition  of  an 
equation  of  state  for  the  high  pressure  and  temperature  regime.  This  allows  the 
model  to  be  used  very  near  an  explosive  source  where  melting  and/or  vaporization 
may  be  important  phenomena.  The  form  of  the  high  pressure  equation  of  state 
adopted  for  the  ARA  model  has  seen  widespread  use  in  the  ground  shock  community. 
Its  development  here  essentially  follows  that  given  by  Shuster  and  Isenberg 
(1972)  and  Schuster  (1981).  At  lower  stresses  (below  the  initiation  of  melting) 
material  behavior  is  controlled  by  the  ARA  three  Invariant  model. 

The  high  pressure  equation  of  state  computes  the  hydrostatic  component  of 
effective  stress,  P,  as  a  function  of  specific  internal  energy,  E,  and  elastic 
volumetric  strain  eve.  Specific  Internal  energy  Is  the  difference  between  heat 
added  to  a  material  and  elastic  volumetric  work  done  by  It,  per  unit  mass.  The 
First  Law  of  Thermodynamics  yields 

dE  =  dQ  -  P  dVe  (19) 

where  dQ  =  Increment  of  heat  added  to  the  substance,  per  unit  mass  and 

dVe  =  elastic  increment  of  specific  volume.  The  volumetric  strain  is  given  by 

the  expression 


V0  -  V  V 

=  Z -  =  1  -  ~  =  1  ~  p0V  (20) 

v0  v0 

so  that 

1 

dVe  - - deve  (21) 

°o 


and  therefore  substitution  of  Equation  (21)  into  Equation  (19)  yields 

P 

dE  =  dQ  +  —  dcve  (22) 

Po 


For  an  Initial  Internal  energy  deposition  at  constant  volume.  Equation  (22) 
yields 

dE  =  dQ  (23) 
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XT<  V  XT*  XT'  V 


and  during  subsequent  adiabatic  deformation.  Equation  (22)  yields 

P 

dE  =  —  dtve  (24 ) 

Po 

The  hydrostatic  stress,  P,  Is  assumed  to  be  the  sum  of  a  solid/liquid 
pressure,  c,  plus  a  vapor  pressure,  p, 

P  =  o  +  p  (25) 

where  the  vapor  pressure,  p,  remains  zero  until  the  specific  internal  energy,  E, 
reaches  the  value,  Em,  required  to  Initiate  melting.  When  E  >  Em,  an  increase 
in  E  contributes  to  the  heat  of  fusion  (melting)  and  vaporization,  as  well  as  to 
Increases  in  both  a  and  p. 

Because  elastic  volumetric  strain,  tve,  and  specific  internal  energy,  E, 
are  the  two  variables  which  determine  the  solid/liquid  pressure,  o,  we  can  write 


The  first  partial  derivative  in  Equation  (26)  is  the  isothermal  elastic  bulk 
modulus,  K, 


The  second  partial  derivative  in  Equation  (26)  is  the  rate  of  increase  of  solid/ 
liquid  pressure  with  respect  to  specific  internal  energy  at  constant  clastic 
volumetric  strain.  It  can  be  expressed  In  terms  of  familiar  quantities  by 
noting  that  the  condition  of  constant  elastic  volumetric  strain  can  be  expressed 
in  the  form 


where  Equation  (27)  yields 


(28) 


(29) 


and  also 


■ 


i 


i 


R 


k :: 

'j 

•< 


k 


R* 


’*  £ 

L.  •* 


V  i 
/ 

«■ 


pj  v?  rj  i 


/3 eve\  /3T\ 

'  »E  K 

'  3T  '3E  , 

(30) 


where  T  *  temperature,  and 


(3£ve\ 

tt) 


■a 


3T  o 
3E> 

3T' 


(31) 

(32) 


The  quantity  <*p  is  the  coefficient  of  thermal  expansion  at  constant  pressure, 
and  the  quantity  Cp  Is  the  specific  heat  at  constant  pressure.  Substitution  of 
Equations  (29)  and  (30)  into  Equation  (28)  yields 


1  ap 

-  da - dE  =  0 

K  Op 


(33) 


so  that 


(-) 

V»E'cve 


K  cq 


(34) 


For  an  initial  energy  deposition  at  constant  volume,  substitution  of  Equation 
(34)  into  Equation  (26)  yields 

K  On  * 


da  = 


dE  =  K  dcve 


(35) 


where 


*  ap 

dcve  =  —  dE 
LP 

During  subsequent  adiabatic  deformation  prior  to  melting,  substitution  of 
Equations  (24),  (25),  (27),  and  (34)  into  Equation  (26)  yields 


(36) 


K  ap  a 

da  K  dcy^  ^  ™  "  dcve  K 


Op  Pq 


[’ 4  &  ’] 


de 


ve 


(37) 


When  E  >  Em,  the  vapor  pressure,  p,  is  calculated  by  an  expression  similar 
to  that  for  adiabatic  compression  of  a  perfect  gas,  for  which 


pV7  =  k 


(38) 


where 

Cp 

r  =  -  (3< 

Cy 

Cy  =  specific  heat  at  constant  volume  and  k  =  constant.  Now  under  adiabatic 
conditions, 


so  that 


dV 

dV 

=  -  p  dV  =  -  pV7  —  = 

-  k  — 

(40) 

V7 

V7 

V 

,  dx 

Vr+11  v  PV 

dE  =  -  k  /  —  =  -  k 

-  =  — 

(41) 

'  X7 

00 

L  -r+1  J  o>  r-1 

p  =  (y-1)  P  E 

(42) 

and  therefore 


The  expression  used  to  calculate  the  vapor  pressure,  p,  is  identical  to  Equation 
(42),  except  a  reduced  specific  energy  is  used  to  account  for  melting  and 
vaporization.  The  equation  is 

p  =  (r-1)  p  E'  (E  >  Em)  (43) 

where 


[  -  ©  1 

E*  =  (E  -  Em)  1  -  e  m  (E  >  Em) 


and  y-1  is  given  by  the  dimensionless  empirical  expression 


7  -  1  =  0.4  +  0.052  In  G  +  0.023  In2  ^ 


where  G  =  p/pQ,  H  =  E  /EQ,  pQ  =  reference  mass  density  (1.0  g/cm  ) ,  and  EQ  = 
reference  specific  internal  energy  (21.171  Te/g).  E*  is  defined  as 


E*  =  Er 


(E*  S  Em) 


E*  =  E1  (E‘  >  Em) 


(46a) 

(46b) 


Equation  (46)  is  necessary  to  avoid  a  logarithmic  singularity  in  Equation  (45). 
It  can  be  shown,  starting  from  Equation  (44),  that  E'  >  Em  when  E  >  2.35  Em. 

The  deviator  stresses  are  reduced  by  the  factor  1  -  E/Em  when  E  <  Em,  so 
when  E  2  Em,  the  material  is  a  fluid  or  a  gas,  and  there  are  no  plastic  strain 
increments.  Thus,  there  is  no  distinction  between  total  and  elastic  volume 
change  whenever  the  vapor  pressure  Is  calculated. 

The  high  pressure  equation  of  state  has  been  implemented  and  Is  included  in 
the  model  version  listed  in  Appendix  A.  The  hydrostatic  behavior  of  the  model 
to  5x10^  MPa  (5  Mbar)  is  shown  in  Figure  5,  both  with  and  without  activation  of 
the  high  pressure  equation  of  state.  Uniaxial  strain  compression  behavior  of 
the  model  with  the  high  pressure  equation  of  state  is  shown  in  Figure  6.  The 
material  melts  at  a  pressure  of  about  2x10$  MPa.  At  this  point  the  deviator 
stresses  (and  the  expansive  yield  surface)  have  been  fully  reduced  to  zero 
(Figure  6b). 
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3.0  COMPUTATIONAL  ISSUES  IN  CONSTITUTIVE  MODELING 


3.1  Tlmestec 


The  timestep  in  time-marching  finite  difference  or  finite  element  solutions 
is  often  determined  based  on  the  Courant  condition.  Simply  put,  this  condition 
permits  a  stress  wave  to  travel  no  more  than  one  zone  thickness  in  one  timestep. 
Thus,  for  a  given  zone: 


At  <  F 


where  At  =  calculational  timestep,  Cp  =  current  compressional  wavespeed, 

Dmin  =  minimum  distance  across  the  zone,  and  F  =  a  safety  factor  (<  1.0)  which 
can  further  restrict  the  timestep.  When  the  Courant  condition  is  employed,  the 
material  model  Is  required  to  report  the  current  wavespeed  In  each  zone.  For  an 
elastic  compressional  wave. 


where  M  =  constrained  modulus  and  p  =  mass  density.  Because  a  zone  may  be 
deforming  plastically,  the  above  elastic  relationship  will  not  always  yield  the 
fastest  signal  propagation  speed.  Therefore,  the  approximate  wavespeed  across  a 
zone  is  then  taken  to  be: 


C-J1 


where  Me  =  the  maximum  of  [Cep(l,l),  Cep( 2,2),  Cep(3,3)]  and  Cep  is  the  6x6 
incremental  elastoplastic  stiffness  matrix.  Using  this  wavespeed  allows  the 
Courant  condition  to  be  closely  followed.  An  additional  safety  factor  of 
F  =  0.9  is  typically  applied. 
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3.2  Numerical  Errors 


The  types  of  numerical  errors  discussed  here  are  primarily  those  which 
occur  within  the  constitutive  model  itself.  Of  the  other  errors  in  a  finite 
difference  or  finite  element  calculation  which  occur  outside  the  material  model, 
those  which  are  controlled  by  artificial  viscosity  are  most  pertinent. 

Artificial  viscosity  is  Intended  to  damp  high  frequency  oscillations.  The 
damping  is  Ignored  when  computing  stress  Increments  within  the  constitutive 
model.  In  conjunction  with  a  rate-dependent  model,  its  use  may  hamper  the 
evaluation  of  strain-rate  effects.  Because  the  ARA  three  invariant  model  is  not 
strain-rate  dependent,  typical  values  of  artificial  viscosity  may  be  employed  to 

smear  shock  fronts  and  limit  grid  oscillations. 

« 

Numerical  errors  produced  in  the  ARA  model  are  primarily  the  result  of  its 
incremental  stiffness  formulation.  The  first  kind  of  numerical  error  is  the 
tendency  for  the  stress  point  to  overshoot  the  expansive  yield  surface  at  low 
values  of  expansive  plastic  work.  Since  this  occurs  when  the  stress  point  is  on 
the  yield  surface  and  pushing  it  out,  it  is  called  "plastic"  overshoot.  Plastic 
overshoot  results  in  a  violation  of  the  expansive  consistency  condition,  which 
states  that 

f '  =  f"  (50) 

P  P 

throughout  yielding.  Plastic  overshoot  is  most  likely  to  occur  as  the  stress 
point  leaves  the  hydrostatic  axis,  because  the  derivative  of  the  expansive 
hardening  function  with  respect  to  expansive  plastic  work  is  infinite  when 

Wp  =  0: 


Because  the  value  of  f^  at  the  end  of  an  Increment  is  computed  from  the  slope 
(dfp/dWp)  at  the  beginning  of  the  Increment,  It  tends  to  be  over-estimated  and 
the  consistency  condition  therefore  violated.  This  phenomenon  can  be  held  In 
check  If  the  strain  Increments  are  kept  very  small  in  this  region.  A  strain 
subcycling  scheme  was  devised  which  evaluates  the  change  in  df"/dW  over  a 


strain  Increment  found  to  have  violated  the  consistency  condition  along  the 
expansive  yield  surface.  This  change  in  slope  is  used  to  break  down  the  total 
current  strain  Increment  into  n  equal  subincrements,  where 


The  model  is  then  internally  cycled  using  these  smaller  Increments.  Compatabil- 
Ity  Is  more  closely  enforced,  eliminating  plastic  overshoot. 

Another  kind  of  error  which  occurs  in  the  ARA  model  (as  well  as  in  other 
models)  Is  associated  with  the  large  strain  increments  which  can  occur  in  one 
timestep  under  explosive  loading.  When  the  stress  point  is  initially  in  the 
elastic  region  of  the  model,  a  large  strain  Increment  can  drive  the  stress  point 
past  one  or  both  yield  surfaces.  There  are  several  ways  of  correcting  this, 
including:  (1)  pulling  the  stress  point  back  to  the  yield  surface,  either  at 
constant  octahedral  normal  stress  or  normal  to  the  yield  surface,  (2)  correcting 
back  at  either  constant  f’c  or  constant  fp,  depending  upon  which  surface  has  been 
violated,  (3)  breaking  down  the  total  strain  increment  into  two  parts  (the  first 
elastic  and  just  sufficient  to  initiate  yielding)  based  on  the  ratio  of  f'/f1', 
or  (4)  breaking  down  the  total  strain  Increment  Into  many  smaller  Increments  and 
subcycling  within  the  model.  This  last  method  Is  currently  employed  in  the  ARA 
model.  Proper  treatment  of  “elastic"  overshoot  was  observed  to  be  very  impor¬ 
tant  in  the  wave  propagation  calculations,  because  of  the  tendency  of  the  stress 
point  to  violate  the  yield  surface  upon  unloading  and  reloading,  especially  from 
a  spalled  condition. 

A  different  kind  of  numerical  problem,  but  not  an  error,  was  encountered 
while  using  the  ARA  model  on  a  VAX  11/750  computer,  which  is  a  much  smaller 
machine  than  the  CRAY,  on  which  the  ARA  model  was  developed.  The  problem  was 
the  occurrence  of  very  large  numbers  associated  with  the  compressive  yield 
surface;  numbers  larger  than  the  computer  could  handle.  Although  the  expansive 
yield  surface  expression, 

/  "'oct  \  /  Pd 

f‘  =  I  -  I  (1  -  f  cos  3u,)  I  -  +  m 

P  V  pa  /  \  °oct 
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has  been  made  dimensionless  through  the  use  of  atmospheric  pressure,  the 
compressive  yield  surface  expression, 

fc  =  3  °oct2  +  rZ  'toct2  (54) 

has  not  been.  So  when,  for  example,  units  of  Pascals  are  being  used  with  stress 
levels  typical  of  blast  loading,  a  quantity  using  (f^)2  will  be  extremely 
large.  The  problem  has  been  circumvented  by  using  units  of  MPa  instead  of  Pa  on 
the  VAX.  There  are  several  ways  to  make  the  compressive  yield  function  dimen¬ 
sionless.  One  possibility  is  to  simply  use  atmospheric  pressure  again  to  cancel 
units: 

3.3  Efficiency 

Evaluating  the  efficiency  of  a  constitutive  model  Involves  answering  three 
questions: 

1.  How  long  does  it  take  the  computer  to  execute  the  model? 

2.  How  much  Information  does  the  model  require  to  be  stored  for  each 
element? 

3.  Are  the  increases  in  run  time  and  storage  required  by  a  more  compli¬ 
cated  model  over  a  simpler  model  offset  by  improved  calculational 
resul ts? 

The  ARA  model  is  a  fairly  long  model  in  terms  of  coding,  as  can  be  seen  in 
Appendix  A.  The  time  it  takes  a  typical  mini -computer  to  run  throuyu  the  uiodei 
Is  compared  with  several  other  models  in  Table  2.  It  is  to  be  expected  that 
more  calculational  steps  will  require  a  somewhat  longer  execution  time.  So,  the 
results  In  Table  2  are  not  surprising.  However,  as  will  be  seen  for  the 
one-dimensional  wave  propagation  calculations,  the  real  run  time  differences 
arise  from  the  strain  subcycling  scheme  used  to  minimize  the  numerical  errors 
discussed  above.  Therefore,  It  Is  not  necessarily  true  that  a  complicated  model 
will  be  significantly  more  expensive  to  run.  What  Is  needed  is  a  more  efficient 
computational  strategy. 
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TABLE  2.  MODEL  EXECUTION  TIMES-SINGLE  ELEMENT  STUDIES 


I 


CPU  Time1  to  Simulate  Laboratory  Test  (sec) 

Model 

Uniaxial 

Strain 

Compression2 

Standard 

Triaxial 

Compression3 

Arbi trary 

Strain 

Path'* 

ARA5 

239 

(ea  =  13.5%) 

301 

34 

CAP 

32 

(ea  =  17.5%) 

48 

6 

AFWL  Engineering 

(ea  =]f7.4%) 

16 

4 

Elastic 

18 

(ea  =  57.0%) 

12 

L 

1  CPU  times  are  to  the  nearest  second  for  the  entire  calculation  with  no 
plotting.  The  computer  used  was  a  Digital  VAX  11/750. 

2  Sample  loaded  to  an  axial  stress  of  40  MPa  and  unloaded  to  an  axial  stress  of 
10  MPa  using  equal  axial  strain  increments  of  0.0025'*.  Note  that  each  model 
resulted  in  a  different  maximum  ea,  as  noted. 

3  Sample  initially  confined  to  3.45  MPa,  loaded  to  102  axial  strain  in  equal 
increments  of  0.00252.  Unloaded  to  0.25  MPa  stress  difference. 

'  Strain  path  shown  in  Figure  35,  simulating  insitu  spherical  wave  propagation. 

‘  Note  that  the  execution  times  achieved  for  the  ARA  model  are  heavily 
influenced  by  the  numerical  correction  scheme  and  can  be  substantially 
improved  by  using  improved  computational  strategies. 


As  coded  in  the  Soil  Element  Model  (SEM),  the  ARA  model  (without  the  high 
pressure  equation  of  state)  has  six  state  variables:  (1)  maximum  past  octa¬ 
hedral  normal  stress,  (2)  compressive  plastic  work,  (3)  expansive  plastic  work, 
(4)  expansive  yield  surface  activation  switch,  (5)  volume  expansion  since  last 
tensile  failure,  and  (6)  Initial  confining  pressure.  This  is  compared  with  no 
state  variables  for  the  elastic  model,  three  for  the  modified  AFWL  Engineering 
model,  and  one  for  the  cap  model.  Each  state  variable  must  be  stored  for  each 
calculatlonal  zone.  The  use  of  six  state  variables  in  the  ARA  model  has  not  yet 
caused  any  storage-related  problems,  either  on  a  large  computer  (CRAY)  or  a 
mini-computer  (VAX  11/750).  Calculatlonal  costs  may  be  increased  due  to  expanded 
core  space  requirements,  but  this  will  not  be  the  case  if  space  is  being 
utilized  which  was  already  available  in  the  code  (as  was  the  case  for  the 
STEALTH  Implementation). 

3.4  Uniqueness  and  Work-Softening 

Uniqueness  In  the  context  of  constitutive  modeling  Is  concerned  with  the 
possibility  of  more  than  one  solution  for  a  given  set  of  stress  or  strain 
conditions.  For  example,  is  It  possible  that  a  total  strain  increment  can 
produce  more  tnan  one  stress  Increment?  If  the  answer  for  a  particular  model  is 
yes,  then  uniqueness  is  violated  and  confidence  in  the  results  generated  by  that 
model  is  greatly  diminished.  For  models  with  two  yield  surfaces  meeting  at.  a 
corner,  the  question  of  uniqueness  at  that  corner  is  particularly  relevant.  The 
ARA  model  employs  a  method  of  choosing  yield  modes  which  was  formulated  to 
insure  uniqueness  [see  Merkle  and  Dass  (1985:  Appendix  I)].  If  a  nonunique 
situation  Is  possible,  the  calculation  Is  stopped.  In  this  way,  a  unique 
solution  has  been  assured  for  all  loading  cases. 

The  tendency  of  some  geologic  materials  to  work-soften,  i.e.,  to  display  a 
decreasing  load  capacity  with  increasing  strain,  has  been  demonstrated  many 
times  in  laboratory  tests.  There  remains  debate  over  the  Interpretation  of 
these  tests  and  whether  or  not  a  soil  sample  is  undergoing  homogeneous  deforma¬ 
tion  at  later  times  In  these  tests.  What  Is  sometimes  Interpreted  as  work¬ 
softening  may  actually  be  a  consequence  of  testing  method,  boundary  conditions, 
or  localized  shear  failure.  However,  there  are  clearly  some  cases  where  soil 
materials  exhibit  a  peak  shearing  resistance  followed  by  a  lesser  residual 
resistance.  The  transition  between  the  two  is  referred  to  as  work-softening. 
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The  expansive  yield  surface  in  the  ARA  model  has  been  formulated  to  account 
for  work-softening.  Currently,  a  material  is  allowed  to  soften  until  zero 
shearing  resistance  remains  (see  Figure  7).  This  is  not  a  good  representation 
for  most  geologic  materials  well  beyond  peak  stress  because  they  tend  to  display 
substantial  residual  shear  strength.  An  expansive  yield  function  which  allows 
residual  strength  has  been  formulated  and  tested  in  the  model  but  is  not  fully 
Implemented. 

During  the  early  stages  of  the  wave  propagation  calculations.  It  was 
observed  that  the  ARA  model  frequently  shut  itself  off  because  It  had  encounter¬ 
ed  the  possibility  of  a  nonunique  solution  in  the  mode  decision  algorithm  at  the 
corner  of  the  yield  surfaces.  This  problem  often  coincided  with  the  onset  of 
work-softening. 

To  expedite  the  continuum  code  checkout  of  the  model,  a  modification  was 
made  to  allow  work-softening  to  be  deactivated.  Thus,  the  yield  surface  can  now 
achieve  a  maximum  value  at  which  it  becomes  stationary  (Figure  7).  The  conse¬ 
quences  of  work-softening  in  dynamic  wave  propagation  certainly  deserve  further 
study,  but  It  was  felt  that  a  complete  treatment  of  this  issue  was  beyond  the 
scope  of  this  effort. 

3.5  Rezoning 

Rezoning  is  an  Important  computational  issue  which  is  often  encountered  in 
finite  difference  calculations  of  blast  and  shock  events  which  employ  a 
tagrangian  grid.  The  need  for  rezoning  arises  when  distortion  of  the  grid 
around  an  explosive  source,  in  a  crater,  or  at  a  stress  concentration  point 
becomes  so  severe  that  It  either  (1)  drives  the  timestep  to  zero,  or  (2)  turns 
the  grid  Inside-out.  Rezoning  is  the  process  of  rearranging  the  grid  (at  one 
instant  of  time  or  over  several  cycles)  so  that  it  Is  again  fairly  uniformly 
spaced,  and  all  zones  resemble  quadrilaterals.  Because  the  numerical  grid  is 
being  remapped  onto  the  material,  as  It  would  be  in  an  Eulerian  grid,  a  calcula¬ 
tion  employing  this  process  is  sometimes  called  arbitrary  Eulerian-Lagranglan 
(ALE).  Rezoning  Is  an  approximate  process  and  does  not  always  conserve  both 
mass  and  momentum.  The  reason  it  involves  (and  Is  dependent  on)  the  constitu¬ 
tive  model  Is  that  as  material  is  transported  and  mixed  from  one  zone  to 
another,  the  state  variables  which  define  each  zone's  material  must  also  be 
consistently  redefined. 
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Figure  8  shows  what  physically  happens  for  one  type  of  Internal  point 
rezone,  using  STEALTH  numbering  conventions.  The  interior  point  common  to  four 
surrounding  zones  is  adjusted  to  a  more  central  location.  In  the  process,  mass 
and  volume  are  exchanged  among  old  zones  to  create  new  zones.  The  quantity  V-jj 
Is  Introduced,  defined  to  be  the  volume  contributed  from  old  zone  1  to  new  zone 
j.  If,  In  creating  new  zone  j,  no  material  is  gained  from  old  zone  1,  then 
=0.  A  matrix  of  weighting  factors  Is  created  for  the  four  new  zones 
created  by  relocation  of  an  Interior  point: 
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F-jj  =  max 
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Fij  Is  defined  as  the  positive  volume  fraction  of  new  zone  j,  which  came  from 
old  zone,  i.  The  "max"  operation  in  Equation  (57)  is  needed  because  V^j  can  be 
negative  due  to  extreme  distortion  before  rezoning.  For  each  rezone  case,  there 
will  be  a  total  of  nine  non-zero  weighting  factors.  Zone  centered  variables, 
such  as  mass  and  internal  energy,  are  then  redistributed  using  these  weighting 
factors.  Material  model  state  variables  must  also  be  distributed  and  a  new 
stress  state  determined  for  each  new  zone.  For  the  ARA  model,  the  parameters 
which  must  be  distributed  are  the  initial  confining  pressure,  the  maximum  past 
pressure,  compressive  and  expansive  plastic  work,  and  volume  expansion  since 
last  tensile  failure  (if  any).  For  adjacent  zones,  redistributing  initial 
pressure  and  reinitializing  confining  pressure  -dependent  properties  Is  not 
critical.  However,  the  new  stress  state  and  state  variables  must  be  consistent 
with  the  new  strain  state,  which  is  a  known  quantity  upon  reconfiguring  the 
zone. 

A  scheme  for  achieving  a  consistent  state  has  not  yet  been  formulated  for 
the  ARA  model,  but  will  be  necessary  for  its  eventual  use  in  ground  shock 
calculations  Involving  severe  environments. 


3.6  Strain  Conventions 


A  large  class  of  finite  difference  codes  (including  TOODY,  HEMP,  STEALTH, 
and  SNEAKY)  do  not  track  displacement  of  grid  points.  Instead,  they  use  current 
strain  rates 
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which,  when  multiplied  by  the  current  timestep  yields 

At 

At  =  —  (59) 

t 

where  t  =  current  zone  length.  These  strains  are  commonly  known  as  true  strains 
and  are  based  on  current  dimensions.  The  properties  of  the  constitutive  models, 
however,  are  typically  formulated  from  laboratory  data  using  engineering 
strains, 

At 

At  =  —  (60) 

*o 

where  to  =  original  zone  length.  The  difference  between  the  two  strains  is 
small  at  small  strains.  But  at  strains  greater  than  about  ten  percent,  the 
difference  becomes  significant  and  can  affect  calculated  results.  Therefore,  it 
has  been  necessary  to  add  variables  in  the  finite  difference  codes  employed  in 
this  study  to  track  original  zone  dimensions  and  to  calculate  actual  engineering 
strains.  These  strains  are  then  used  in  the  material  models  anc  are  consistent 
with  the  development  of  the  models  and  their  parameters. 


4.0  WAVE  PROPAGATION  CALCULATIONS 


4.1  ARA  Model  Implementation 

The  ARA  model  has  been  implemented  in  a  fashion  compatible  with  finite 
difference  or  finite  element  code  applications.  The  model  Is  formulated  to 
operate  under  strain  control,  where  the  total  strain  increment  is  known  and  the 
resultant  stress  Increment  Is  determined  by  the  model.  The  incremental  elasto- 
plastic  stiffness  approach  used  by  the  ARA  model  [Merkle  and  Dass  (1985)],  is 
fundamentally  different  than  the  trial  and  error  failure  surface  correction 
procedure  employed  by  many  current  models,  including  the  AFWL  Engineering  and 
cap  models.  The  Incremental  stiffness  procedure  more  accurately  tracks  plastic 
strains  and  plastic  work,  although  some  errors  are  produced  by  extrapolating 
stiffness  from  an  old  stress  state  to  a  new  one  (see  Section  3.2). 

Because  the  ARA  model  was  developed  using  the  Soil  Element  Model  (Dass, 
Bratton,  and  Higgins  (1981)],  its  Implementation  has  kept  pace  with  its  improve¬ 
ment  and  change.  The  model  has  been  extensively  tested  in  a  single  element  mode 
under  many  kinds  of  laboratory  stress  and  strain  paths.  Applying  the  model  to 
wave  propagation  problems,  however,  did  require  some  modification,  lhe  behav¬ 
ioral  improvements  are  discussed  in  Section  2,  ana  numerical  improvements  are 
discussed  in  Section  3. 

4.2  Initial  Anisotropic  Stress  State 

The  stress  field  at  a  point  in  an  earth  mass  initially  at  rest  under  the 
force  of  gravity  depends  on  depth,  local  tectonic  conditions,  and  material 
properties.  If  the  action  of  gravity  during  geologic  history  is  idealized  as  a 
uniaxial  strain  compression  process,  then  calculating  Lite  insitu  stress  field 
may  precede  accordingly.  The  calculation  is  relatively  easy  for  an  elastic 
material . 

The  elastic  geostatic  octahedral  normal  stress  is 


where  oz  is  the  vertical  effective  stress  at  the  depth  of  interest  due  to 
overburden,  K0  is  the  coefficient  of  earth  pressure  at  rest,  and  v  is  Poisson's 
ratio.  The  insitu  stress  state  is  then  a  point  on  the  elastic  uniaxial  stress 
path  shown  in  Figure  9a. 

For  an  elastic-plastic  material,  which  can  actually  be  at  incipient  shear 
failure  under  insitu  stress,  or  for  which  there  can  be  model  state  parameters 
which  cannot  be  directly  determined,  it  is  possible  to  initially  load  each 
element  In  uniaxial  strain  compression  (Figure  9b).  At  the  proper  vertical 
stress,  each  element  is  then  at  equilibrium  under  gravity  with  correct  hori¬ 
zontal  stresses  as  well  as  state  parameters.  If  one  or  more  model  parameters 
depend  on  "initial"  confining  pressure,  however,  this  process  becomes  more 
complicated.  These  model  parameters  are  used  to  load  the  element  uniaxially, 
but  the  final  model  parameters  actually  depend  on  the  at-rest  anisotropic 
stress  state.  If  these  parameters  cannot  be  explicitly  determined,  it  would 
appear  that  some  form  of  iteration  is  necessary  to  arrive  at  the  true  "initial" 
condition. 

The  ARA  model  has  several  parameters  which  depend  on  initial  confining 
pressure.  These  parameters  control  the  work  hardening  functions  and  the 
expansive  plastic  potential  surface,  as  well  as  the  initial  elastic  moduli.  In 
order  to  characterize  these  parameters  correctly,  a  procedure  has  been  devised 
to  approximate  initial  anisotropic  consolidation.  This  procedure  is  exercised 
for  each  depth,  so  that  as  depth  increases,  each  zone  will  have  unique  initial 
conditions.  The  steps,  shown  in  Figure  9c,  are  as  follows: 

(1)  Eliminate  the  usual  model  initialization  which  occurs  when  parameters 
are  input.  Or  initialize  the  model  parameters  to  a  very  low  isotropic 
pressure,  approximately  one-tenth  of  an  atmosphere  (point  a  in 
Figure  9c). 

(2)  Calculate  the  stress  state  which  would  be  achieved  due  to  overburden 
at  the  depth  of  interest  if  the  material  were  elastic  (as  was  dis¬ 
cussed  above).  Use  the  unloading-reloading  Poisson's  ratio,  vu,  to  do 
this  (from  Point  a  to  Point  b). 

(3)  Reinitialize  the  model  parameters  using  the  octahedral  normal  stress 
found  in  Step  2  (point  c). 
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Figure  9.  Application  of  Geostatic  Stresses. 


(4)  Load  the  model  In  uniaxial  strain  compression  until  the  vertical 
stress  reaches  the  overburden  level  (from  point  c  to  d).  Save  the 
Insltu  stress  state  and  the  value  of  the  expansive  plastic  work,  Wp, 
at  this  point. 

(5)  Reinitialize  all  other  model  parameters  using  the  pressure  level 
reached  in  Step  4  (point  e). 

Achieving  an  appropriate  and  consistent  insltu  stress  state  is  Important  to 
a  successful  calculation  for  several  reasons.  First,  the  calculational  grid 
will  be  stable  under  gravity  forces  prior  to  arrival  of  a  stress  wave.  Second, 
the  material  behavior  of  models  such  as  the  ARA  conic  model  can  be  quite 
sensitive  to  initial  confining  pressure.  And  third,  the  orientation  of  the 
stress  wave  with  respect  to  the  orientation  of  the  insitu  stresses  is  important 
when  the  dynamic  and  static  stresses  are  of  the  same  order  of  magnitude.  Figure 
10  illustrates  the  second  and  third  of  these  effects.  Shown  are  the  stress 
paths  due  to  a  hypothetical  insitu  spherical  strain  path.  The  strain  path  used 
here  is  the  same  as  that  used  before  in  Figure  3.  Three  cases  were  calculated: 

(1)  Initial  Isotropic  compression  to  6  MPa.  Since  the  initial  stresses 
are  equal  in  the  x,  y,  and  z  directions,  direction  of  load  application 
does  not  affect  this  case. 

(2)  Initial  anisotropic  consolidation  to  ox  =  6  MPa.  T he  strain  path  is 
applied  as  if  the  stress  wave  were  traveling  in  the  vertical  (x) 
direction.  Behavior  is  qualitatively  similar  to  Case  (1)  but  quanti¬ 
tatively  quite  different. 

(3)  Initial  anisotropic  consolidation  to  ox  =  6  MPa  (same  as  Case  (2)). 

The  strain  path  is  applied  as  if  the  stress  wave  were  traveling  the 
horizontal  (y)  direction.  Note  the  difference  in  initial  behavior  due 
to  the  initial  drop  in  shear  stress.  This  same  response  will  be 
apparent  in  Section  4.5  for  the  two-dimensional  DIHEST  calculation. 


4.3  Description  of  Modeled  Soil 

The  ARA  model  is  intended  to  be  quite  general  and  is  capable  of  modeling 
many  types  of  soil,  as  well  as  other  kinds  of  materials.  Fundamental ly  differ¬ 
ent  types  of  soil  response  can  be  matched  through  the  parameter  determination 
process.  The  soil  modeled  in  the  wave  propagation  problems  discussed  here  is  a 
dry  alluvium  representative  of  the  alluvial  materials  found  across  the  desert 
southwest  and  the  basin  and  range  topographies  of  the  United  States.  Dry 
alluvium  was  chosen  for  several  reasons.  First,  dry  alluvium  is  of  great 
interest  to  the  Air  Force  because  of  the  ISST  test  program  currently  being 
conducted  in  Yuma,  Arizona.  Secondly,  there  is  a  good  deal  of  data  (both 
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(a)  Isotropic  consolidation  to  P  =  6  iPa 

(b)  Anisotropic  to  J  =6  MPa  -  ' 
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a.  Stress  Path  Response. 
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b.  Pressure-Volume  Response. 

Figure  10.  Some  Effects  of  Anisotropic  Consolidation  in  the  ARA  Model 
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laboratory  and  Insitu)  available  from  many  test  programs  performed  in  similar 
materials.  Third,  a  dry  material  was  necessary  to  avoid  the  additional  computa¬ 
tional  issues  Involved  when  propagating  waves  through  saturated  or  partially 
saturated  media.  Finally,  the  dry  alluvium  used  here  Is  the  same  material  used 
to  demonstrate  the  model  against  laboratory  data  by  Merkle  and  Dass  (1985). 

Table  3  shows  the  ARA,  CAP,  AFWL  Engineering,  and  elastic  model  prameters 
used  for  the  calculations  reported  here.  The  fact  that  data  from  remolded 
laboratory  samples  was  used  to  fit  the  model  parameters  Is  important  because  It 
limits  the  potential  accuracy  of  the  model  in  predicting  Insitu  test  results. 
Even  when  undisturbed  samples  are  used  to  fit  a  model,  the  predicted  response 
using  a  laboratory  model  is  commonly  different  from  the  dynamic  insitu  response, 
for  both  simple  and  complex  geometries.  When  simple  models  are  used,  the 
parameters  can  be  adjusted  to  yield  a  best  estimate  of  what  Insitu  response  will 
be.  Figure  11  shows  preliminary  estimates  of  laboratory  and  insitu  uniaxial 
strain  compressibility  for  ISST  alluvium  [Jackson  (1984)].  No  adjustment  of 
this  kind  has  been  utilized  for  this  study.  The  result  will  be  poor  agreement 
between  the  calculations  and  the  data  from  insitu  events.  Because  the  observed 
insitu  responses  are  typically  stiffer  (at  low  to  Intermediate  stress  levels) 
than  the  laboratory  response,  calculated  peak  motions  based  on  laboratory- 
derived  properties  will  be  too  high  and  peak  stresses  will  be  too  low. 

Note  also  that  although  only  one  material  was  used  for  all  the  wave 
propagation  calculations,  the  materials  in  which  the  tests  were  conducted  varied 
significantly.  All  the  Insitu  tests  were  in  dry  alluvium  but  there  were 
variations  from  site  to  site,  and  with  depth  within  each  site.  Because  the 
purpose  of  these  calculations  was  to  check  out  the  ARA  model  in  continuum  code 
problems,  however,  one  material  was  used  for  all  tests  at  all  depths.  This 
allowed  a  direct  comparison  of  model  responses  for  different  geometries  ana 
insitu  conditions. 


4.4  One-Dimensional  Calculations 


4.4.1  Code  Description.  The  finite  difference  code  used  for  the  one- 
dimensional  calculations  Is  an  adaptation  of  the  mechanical-only  portion  of 
SNEAKY,  written  by  Hart  (1981).  It  has  been  Incorporated  In  the  ARA  Soil 
Element  Model  as  a  boundary  condition  option  along  with  the  laboratory  test  and 
arbitrary  strain  path  options.  The  frame  of  reference  in  the  code  is 
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TABLE  3.  REMOLDED  CARES-DRY  ALLUVIUM  CONSTITUTIVE  MODEL  PROPERTIES 


ft* 

SjmPOl 

Parameter 

iLM  home 

Value 

Uni  l  v 

■ 

(1)  AAA  Model 

1 

Hur 

Clastic  Modulus  Coefficient 

AX  UN 

363.  5 

_ 

5\ 

R 

Clastic  Modulus  Exponent 

AN 

0.8412 

- 

VO 

Clastic  Poisson's  Ratio 

APOI 

0.20 

- 

- 

Unload-Reload  Hysteresis  Switch 

AHSUITCH 

0.0 

- 

54 

- 

No.  Collapse  Function  Segments 

ACRV 

3.0 

- 

a 

- 

Work  Softening  Switch 

AWSOFT 

0.0 

- 

f*Z 

Cl 

AACC(l) 

4.645E-5 

- 

pi 

AAPC(l) 

1.401 

- 

C2 

Compressive  Hardening  Constant 

AACC ( 2 ) 

6.086C-2 

- 

m 

P2 

Compressive  Hardening  exponent 

AAPC ( 2 ) 

0.4667 

- 

y 

C3 

AACC ( 3 ) 

0.7616 

- 

P3 

AAPC ( 3 ) 

0.2688 

- 

r 

elliptical  Cap  Shape 

AR 

0.25 

- 

£ 

Expansive  Yield  Constant 

ACY 

0.1111 

- 

v 

• 

expansive  Yield  Constant 

AMY 

2.876C-4 

- 

> 

el 

expensive  Failure  Constant 

ACTA1 

0.6454 

- 

P 

expansive  Hardening  Constant 

APBAR 

0.5067 

- 

% 

CxpansIve  Hardening  exponent 

AL 

0.8691 

- 

• 

Cxpansive  Hardening  Constant 

AALPH 

5.000 

- 

( 

CxpansIve  Hardening  Constant 

ABCTA 

-2.631E-3 

- 

t 

Plastic  Potential  Constant 

ATG 

-0.9646 

- 

R 

Plastic  Potential  Constant 

ARG 

2.182E-3 

- 

■ « 

s 

Plastic  Potential  Constant 

ASG 

1.B60 

- 

T 

Tensile  Strength 

APCX 

0.0 

Pa 

*  * 

L  • 

9 

Mass  Density 

RHORCF 

1900. 

kg/m3 

(2)  Cap  Model 

XI  1 

AXI 

4.0E9 

XI  } 

Bulk  Modulus  Parameters 

AK1 

0.0 

X2  ) 

Ax  2 

0.0 

61  1 

AG! 

3.0E9 

61  i 

Shear  Modulus  Parameters 

AG1 

0.0 

62  ) 

AX  2 

0.0 

C  ) 

AC 

0.288E6 

" 

Failure  Surface 

AM 

CCC 

0.215 

0.0 

6  ) 

88 

0.0 

R1  1 

ARI 

2.5 

R1  | 

Cap  Shape 

AR1 

0.0 

R2  ) 

AR2 

0.0 

s  I 

Cap  Hardening  Parameters 

AH 

A0 

0.200 

0.018E-6 

9 

Mass  Density 

RH0REF 

1900. 

(3) 

AFHL  Engineering  Model 

a 

Mass  Density 

RH0REF 

1900. 

kg/m3 

T 

Tension  Cutoff 

ST1 

-0.288E6 

Pa 

Y 

Yield  Intercept 

Y1 

0.288E6 

Pa 

S 

Failure  Surface  Slope 

SI 

0.215 

- 

VM 

von  Mlses  Cutoff 

VM1 

175E6 

Pa 

Hydrostat  (No.  Load  Slopes  •  8,  No. 

Unload  Slopes  •  5) 

Loading  Scymnts: 

1 

2 

3  4 

5 

6 

7 

8 

XA 

Bulk  Modulus  (Pa) 

9.302E7 

6.261E7 

1.491E8  3.530E8 

1.088E9 

3.419E9 

9.042E9 

7.000E11 

«k 

Strain  Breakpoint 

0.001181 

0.008191 

0.1292  0.1642 

0.2014 

0.2294 

0.2516 

1.0000 

*> 

Poisson's  Ratio 

0.32 

0.32 

0.32  0.32 

0.32 

0.32 

0.32 

0.32 

Unloading  Salients: 

1 

2 

3  4 

5 

Hu 

Bulk  Modulus  (Pa) 

9.000E11 

4.500E11 

1.39E10  4.725E9 

1.000E9 

pu 

Pressure  Breakpoint 

3.7EB 

2.8E8 

3. 0E  7  2.0E7 

-0.288E6 

vu 

Poisson's  Ratio 

0.20 

0.20 

0.20  0.20 

0.20 

f  ■ 

[  o 

(4)  rustic  Model 

K  Bulk  Modulus 

BUix 

B0E6 

Pa 

-  | 

♦ 

1 

* 

6  Shear  Modulus 

SHEAR 

31 E  b 

Pa 

a  Mass  Density 

SHORE f 

1900. 

kg/m3 

35 
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Lagrangian,  so  no  material  is  transported  from  zone  to  zone.  (No  rezoning  was 
used  for  these  calculations. )  The  calculations!  sequence  for  a  time  step  is 
shown  in  Figure  12.  The  overall  formulation  of  SNEAKY  is  very  similar  to 
STEALTH  which  was  used  for  the  two-dimensional  calculations  discussed  in 
Section  4.5 

Only  quadratic  artificial  viscosity  was  used  for  the  one-dimensional 
calculations.  Since  quadratic  artificial  viscosity  was  activated  only  during 
compression,  the  shock  fronts  were  spread  out  but  the  zone-to-zone  numerical 
oscillations  seen  after  the  peak  were  not  damped.  The  artificial  viscous  stress 
in  SNEAKY  is  given  by 


(3x  \  ? 

—  1  (63) 

where  p  =  current  mass  density,  h  =  zone  thickness,  3x/3x  =  strain  rate  in  the 
direction  of  propagation,  and  Cq  =  dimensionless  constant.  Cq  was  set  equal  to 
2.0  for  these  calculations. 

4.4.2  Planar  (HEST)  Calculations.  The  HEST  (High  Explosive  Simulation 
Technique)  is  commonly  used  for  simulating  superseismic  airblast  effects  from  a 
nuclear  detonation.  Within  the  working  volume  and  simulation  time,  which  depend 
on  the  extent  of  the  HEST  cavity,  the  simulator  produces  loading  conditions 
which  are  essentially  one-dimensional  uniaxial  strain  compression.  Figure  13 
shows  the  experimental  layout  for  SIMCAL  3  (Simulation  CALibration  3),  which  was 
a  HEST  test  performed  by  the  Air  Force  at  the  HAVE  HOST  test  site  near  Yuma, 
Arizona  in  1979  [AFWL  (1981)].  This  test  was  chosen  for  analysis  here  because 
it  was  fielded  adjacent  to  the  ISST  test  site,  and  because  it  was  a  fairly 
successful  test  in  terms  of  data  recovery. 

The  calculational  idealization  of  this  test  is  shown  in  Figure  14.  One 
material  was  used  for  the  entire  calculational  grid.  Layering,  which  was 
observed  at  the  site,  has  been  neglected.  The  mesh  consisted  of  99  grid  points 
with  fairly  fine  zoning  near  the  surface.  The  bottom  boundary  was  deep  enough 
to  avoid  reflections  within  the  time  of  interest.  Acceleration  due  to  gravity 
was  Included,  and  the  grid  was  Initially  loaded  to  an  anisotropic  stress  state. 
Target  point  depths  were  chosen  to  correspond  to  data  locations  in  SIMCAL  3. 


For  All  Grid  Points 


1.  Calculate  Acceleration 


x-jn  Using  Conservation  of  Momentum 


2.  Calculate  Velocity 

i^n+1/2  Integration  of  Acceleration 

I 

3.  Calculate  Displacement 

-  xin+l  Integration  of  Velocity 


l 

Calculate  Strain 

En+1  from  Velocity  Gradient 


Calculate  Stress 
°1jn+^  ^rom  SEM  Material  Model 

l 

Calculate  Time  Step 
-  Atn+^  From  Stability  Criterion 


Figure  12.  SNEAKY  1-0  Calculatlonal  Sequence  for  A  Time  Step 

(Hart  (1980:3-27)] 


b.  Elevation. 


Figure  13.  SIMCAL  3  Experiment  Configuration. 
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Figure  14.  1-D  Planar  Wave  Propagation  Cal culational  Configuration. 
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The  exponentially-decaying  pressure  waveform  shown  In  Figure  15  was  applied 
to  the  ground  surface.  It  corresponds  to  a  424  kt,  313  m  range  nuclear  airblast 
loading  [Brode  and  Speicher  (1984)].  This  Is  an  impulse  fit  to  the  recorded 
HEST  overpressures  In  SIMCAL  3  as  developed  by  the  Air  Force  Weapons  Laboratory 
(1981).  The  calculation  was  run  out  to  40  ms,  at  which  point  the  impulse  had 
not  yet  reached  its  peak.  (Impulse  peaks  at  about  1.2  seconds  for  this  yield 
and  range. ) 

Calculated  vertical  velocity  and  displacement  at  five  locations,  using  the 
ARA  model,  are  shown  in  Figures  16a  and  b,  respectively.  It  is  very  clear  that 
in  this  highly  compactive  media,  the  entire  grid  moves  downward  nearly  as  a 
rigid  body  behind  the  wave  peak.  Figures  17a  and  b  show  the  stress  path  and 
pressure-volume  response,  respectively,  for  three  depths.  Because  the  loading 
is  uniaxial  strain  compression,  this  kind  of  response  is  entirely  predictable 
from  the  laboratory  uniaxial  test  simulations  run  previously  with  this  model. 

The  calculated  results  do  not  predict  the  test  results  very  well,  because 
the  parameters  were  not  fit  to  match  the  insitu-based  loading  hydrostat.  This 
is  shown  in  Figure  18,  where  it  is  seen  that  the  calculated  velocity  peaks  are 
too  high,  the  arrivals  too  slow,  the  rise  times  too  fast,  and  the  duration  too 
long. 

The  identical  problem  was  run  using  three  other  models:  the  AFWL 
Engineering  model,  the  cap  model,  and  an  elastic  model.  A  typical  comparison 
between  velocity  waveforms  is  shown  in  Figure  19  for  two  depths.  The  elastic 
model  results  are  fundamental ly  different  because  there  is  no  permanent  com¬ 
paction.  The  remaining  models  all  give  very  similar  results.  The  ARA  model  is 
a  little  softer  than  the  AFWL  or  cap  models,  as  is  seen  in  Figure  20b.  This  is 
partly  due  to  its  initialization  at  a  very  low  confining  pressure  at  this  depth. 
Figure  20a  demonstrates  that  all  models  give  similar  loading  stress  paths.  And, 
since  they  all  have  identical  unload-reload  Poisson's  ratios,  all  three  com¬ 
pactive  models  have  Identical  unloading  stress  paths.  Figure  21  compares  peak 
velocity  attenuation  between  the  models. 

4.4.3  Spherical  (Burled  HE  Sphere)  Calculations.  The  purpose  of  perform¬ 
ing  one-dimensional  calculations  in  geometries  other  than  planar  is  to  more 
fully  exercise  the  ARA  model  prior  to  delving  into  the  two-dimensional  realm  of 
arbitrary  stress  and  strain  paths.  Planar  wave  propagation  subjects  the  model 


Pressure  (MPa) 


ir»  V  »  *7*  L-w  L>  w  i"V  V-  -  W>i  .>  *?%  W  y%'  „"W  VV  J^T ^T\  l VLVv^T\  iVln.Tr w^n.  •V*' w**- 


Brode-Spei< 


W  =  424  kt 

R  =  313  m 

P  =  9.89  MPa 
o 


Time  ( s ) 


Figure  15.  Applied  Surface  Pressure  and  Impulse  Histories 
for  the  1-D  Planar  HEST  Problem. 
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Figure  20.  Model  Behavior  Comparisons  for  1-D 
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Figure  21.  1-D  Planar  Peak  Velocity  Attenuation  -  Model  Comparisons 
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to  strictly  uniaxial  strain  conditions  (two  principal  strains  equal  to  zero). 
Spherical  wave  propagation  is  the  next  step,  where  two  principal  strains  (hoop) 
are  also  equal,  but  not  generally  zero.  The  resultant  stress  paths  are  more 
general,  and  exercise  the  shear  yield  surface  and  tensile  failure  components  of 
a  model  much  more  than  does  uniaxial  strain. 

The  problem  selected  for  these  calculations  was  the  detonation  of  a  small 
fully-buried  sphere  of  high  explosive.  This  kind  of  experiment  has  been 
conducted  many  times  in  conjunction  with  several  DOD  programs;  e.g.,  MOLE, 

ESSEX,  and  most  recently,  ISST.  The  MAT  PROP  #1  charge  was  a  twenty  pound  (9.05 
Kg)  sphere  of  C-4  explosive,  buried  twenty  meters  deep  in  alluvium  at  the  ISST 
test  site  near  Yuma,  Arizona  [Trulio  (1983)].  Figure  22  shows  this  configura¬ 
tion.  The  test  itself  was  not  very  successful  (data  recovery  was  minimal),  but 
it  does  provide  a  case  of  a  small  buried  charge  in  alluvium.  Figure  23  shows 
the  calculational  grid  for  these  calculations.  A  uniform  geology  was  assumed. 

An  initial  Isotropic  prestress  corresponding  to  a  depth  of  20  m  was  applied,  but 
no  gravity  forces  were  used  in  the  calculation.  The  HE  sphere  was  modeled  using 
the  JWL  equation  of  state  [Lee,  et  al.  (1968)]  for  ni tromethane.  The  JWL 
equation  of  state  is  an  empirical  relationship  used  to  predict  the  behavior  of 
explosives  by  accounting  for  large  expansion  of  the  detonation  products.  The 
equation  for  pressure  (P)  is: 


P  =  A 


rRl  V  + 


n  1  / 

,"r2  *  +  _ 

V 


where  V  stands  for'  relative  volume  (V/V0),  E  is  energy  density  (per  unit 
volume),  and  A,  B,  R] ,  R^,  and  w  are  experimental ly  derived  constants.  Table  4 
lists  the  parameters  for  ni tromethane,  which  were  used  to  approximate  C-4.  The 
explosive  was  not  burned  (i.e.,  no  detonation  front  was  propagated),  so  peak 
cavity  stress  was  lower  than  would  actually  occur.  However,  actual  peak  cavity 
stress  decays  so  rapidly  with  range  that  it  is  inmediately  damped  down  to  a 
non-burned  level  in  a  calculation  using  these  zone  sizes  (0.1  m).  Burning  the 
explosive  would  therefore  have  been  Inconsequential.  Figure  24  shows  the 
pressure  history  generated  In  the  cavity. 

The  spherical  calculations  were  run  to  only  a  very  short  time,  about  2  ms. 
This  is  shorter  than  practical  for  normal  production  calculations  of  this  type, 
but  was  desirable  here  for  two  reasons.  First,  results  from  every  cycle  could 
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TABLE  4.  JWL  EQUATION  OF  STATE  PARAMETERS  FOR  NITROMETHANE 
[Lee,  et  al.  (1968:5)] 


Parameter  j 

Symbol 

j  Value  j 

Units 

A 

2.0925x10^1 

J/M3 

8 

5.0895xl09 

J/M3 

Constants 

R1 

4.4 

— 

r2 

1.2 

— 

to 

0.30 

— 

Mass  Density 

Po 

1128. 

kg/m3 

Total  Available 

Energy 

E0 

5.1xl09 

J/M3 
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be  plotted,  which  showed  exactly  how  the  model  was  performing  without  any  gaps 
due  to  plotting  restrictions.  Second,  the  zones  immediately  adjacent  to  the 
source  were  of  most  Interest  because  of  their  severe  compression  and  then  rapid 
hoop  expansion.  It  is  often  true  that  for  very  severe  environment  calculations, 
the  first  few  cycles  tell  the  whole  story.  Running  this  calculation  out  to 
longer  duration  would  also  eventually  require  rezoning,  which  has  not  yet  been 
fully  developed  for  the  ARA  model. 

Radial  velocity  and  displacement  waveforms  generated  using  the  ARA  model 
are  shown  in  Figure  25.  Very  steep  motion  gradients  close  to  the  source  are  due 
in  large  part  to  grid  size  effects.  At  a  loading  wavespeed  of  roughly  400  m/s, 
0.1  m  zones  will  act  as  a  2000  Hz  low  pass  filter  on  the  pressures  generated  in 
the  cavity.  Velocity  waveforms  are  compared  with  results  from  several  other 
models  in  Figure  26.  The  elastic  model  is  the  only  model  which  produces  a 
substantially  different  waveform,  mainly  because  no  shear  failure  occurs.  This 
is  shown  in  Figure  27,  where  stress  paths  at  the  0.44  m  target  point  are 
overlaid  for  the  different  models.  The  highest  shear  stress  is  achieved  during 
the  initial  loading.  Subsequent  loops  in  the  ARA  and  AFWL  Engineering  model 
stress  paths  are  caused  by  rejoin  after  spall.  Spall  is  caused  by  the  rapid 
hoop  expansion,  rejoin  is  caused  by  the  relatively  high  pressure  which  remains 
in  the  cavity.  Figure  28  compares  the  strain  paths  calculated  using  the  four 
model s . 

4.4.4  Cylindrical  (CIST)  Calculations.  Within  the  constraints  of  one¬ 
dimensional  wave  propagation,  cylindrical  geometry  is  capable  of  producing  the 
most  general  stress  paths.  This  is  because  the  three  principal  stresses  and 
strains  are  not  equal.  Thus,  a  material  model's  ability  to  account  for  many 
situations  encountered  in  two-dimensional  calculations  is  tested. 

The  Cylindrical  Insitu  Test  (CIST)  geometry  was  chosen  for  these  calcula¬ 
tions  because  of  the  large  insitu  dynamic  data  base  generated  by  these  experi¬ 
ments.  There  have  been  twenty-three  CIST’s  to  date,  in  dry  soil,  wet  soil,  and 
rock  geologies.  All  have  had  essentially  identical  charge  design:  a  two  foot 
diameter  borehole  filled  with  racked  400-grain  PETN  detonating  cord,  and  an 
explosive  density  of  5  lbs/linear  ft  of  cavity.  The  nominal  peak  cavity 
pressure  for  this  explosive  configuration  is  about  40  MPa.  Subsequent  decay  of 
cavity  pressure  varies  with  depth  and  the  properties  of  the  surrounding  geologic 
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material.  CIST  18  was  conducted  at  the  HAVE  HOST  site,  at  a  location  close  to 
the  two  experiments  discussed  previously.  It  was  actually  two  tests,  CIST  18  s 
(shallow)  and  CIST  18  d  (deep).  The  configuration  of  CIST  18  s&d  is  shown  in 
Figure  29.  More  Information  is  given  by  Amend,  Ullrich,  and  Thomas  (1977). 

The  axisynmetric  calculation  used  to  model  CIST  18  Is  shown  In  Figure  30. 
One  material  was  used  throughout  the  grid.  Gravity  was  not  used,  but  an 
isotropic  prestress  of  0.18  MPa,  corresponding  to  about  a  10  m  depth,  was 
applied.  The  target  points  used  include  the  typical  CIST  gage  radii  of  0.91, 
1.52,  and  2.44  meters. 

The  appropriate  pressure  boundary  for  driving  CIST  calculations  is  always 
uncertain  because  there  have  been  few  successful  measurements  made  in  the 
explosive  cavity.  CIST  18  was  one  of  the  few  tests  which  yielded  a  reasonable 
pressure  history,  so  this  was  used  to  fit  a  two-term  exponential  function.  The 
form  of  the  equation  used  was: 

P(t)  =  21  e-900t  +  19  e"165t  (t  >  0.000125)  (65) 

with  units  of  MPa  and  seconds.  This  pressure  history  and  its  impulse  are  shown 
in  Figure  31.  These  calculations  were  run  to  a  duration  of  about  40  ms. 

Calculated  motions  using  the  ARA  model  are  shown  in  Figure  32.  A  tendency 
for  low  frequency  outward  flow  is  noted  at  all  ranges  due  to  shear  failure  and  a 
very  high  unload-reload  modulus.  The  higher  frequency  motions  superimposed  on. 
this  are  caused  by  post-spall  rejoin  signals  generated  near  the  source.  These 
late  time  spikes  are  not  real,  i.e.,  they  are  not  observed  in  test  data,  and  are 
consequences  of  using  a  minimum  pressure  cutoff  for  a  tensile  failure  criterion. 
Figure  33a  shows  the  volume  compression  response  calculated  at  several  ranges. 
The  tendency  for  the  model  to  continue  to  compress  during  post-peak  cycles  of 
reloading  is  the  most  interesting  aspect  of  this  behavior.  Figure  33b  shows 
strain  paths  for  these  same  ranges. 

Calculated  velocity  is  compared  with  some  representative  composite  data 
from  CIST  18  s&d  in  Figure  34.  Again,  the  lab-based  model  is  clearly  deficient 
in  predicting  insltu  motions. 

Cylindrical  results  were  generated  using  the  three  other  models  for 
comparison  with  the  ARA  model.  The  velocity  waveforms  for  the  ARA  and  the  AFWL 
Engineering  models  are  very  similar  (Figure  35).  The  cap  model  seems  to  do 
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better  and  does  not  tend  to  flow  out  as  do  the  others.  The  reason  Is  its  volume 
response  under  these  conditions,  shown  in  Figure  36.  The  cap  model  has  a  much 
higher  effective  unload-reload  modulus  and  does  not  cycle  as  much  as  either  the 
ARA  or  AFWL  Engineering  models.  Strain  paths  at  two  ranges  are  compared  in 
Figure  37. 

4.4.5  Discussion  of  One-Dimensional  Calculations.  In  an  effort  aimed 
primarily  at  successful  model  implementation,  emphasis  is  naturally  on  computa¬ 
tional  aspects  rather  than  on  actual  soil  behavior.  So  the  first  thing  to  be 
said  about  the  one-dimensional  results  is  that  they  show  the  ARA  model  is 
capable  of  functioning  In  a  dynamic  finite  difference  wave  propagation  environ¬ 
ment.  The  fact  that  the  ARA  model  results  do  not  always  match  observed  insitu 
behavior  is  not  too  troublesome  at  this  point.  The  data  comparisons,  when 
shown,  provide  a  goal  toward  which  the  model  fitting  process  can  precede. 

Most  of  the  effort  which  went  into  the  calculations  with  the  ARA  model  was 
spent  on  eliminating  the  kinds  of  numerical  errors  discussed  in  Section  3.2. 
Yield  surface  violation  due  to  both  plastic  and  elastic  overshoot  was  the  most 
serious  problem.  The  subcycling  technique  was  devised  to  circumvent  reformula¬ 
tion  of  the  expansive  plastic  work  function  to  avoid  the  singularity  at  the 
isotropic  axis,  or  modification  of  the  solution  technique  (from  incremental 
stiffness  to  predictor-corrector) .  Subcyling  may  not  be  the  most  efficient 
approach.  The  next  biggest  problem  turned  out  to  be  the  tensile  failure  loqic 
(spall  model).  With  the  tendency  of  the  stress  point  to  move  auickly  down  the 
yield  surface,  many  the  zones  (particularly  near  the  source)  spent  most  of 
the  time  at  the  apex  of  the  expansive  yield  surface.  Tracking  spall  strain  and 
providing  for  orderly  rejoin  was  therefore  critical  for  overall  success. 

Gravity  initialization  was  the  third  major  outcome  of  the  one-dimensional 
calculations. 

The  ARA  model  was  the  most  time-consuming  of  the  four  models  compared. 

Table  5  compares  times  spent  in  various  parts  of  the  calculations.  A  grid  cycle 
is  the  calculation  of  the  response  for  the  entire  grid  for  one  time  step.  Grid 
cycle  times  varied  for  the  ARA  model  and  the  cap  models  because  a  different 
number  of  zones  may  be  subcycled  for  any  given  cycle.  Since  subcycling  occurs 
mostly  near  the  Isotropic  axis  for  the  ARA  model,  application  of  anisotropic 
gravity  stress  (grid  setup)  in  the  1-D  planar  case  took  a  great  deal  of  time, 
but  substantially  reduced  subsequent  dynamic  calculation  time. 


.vw  .V  -  . 


— —  ARA  Model 

. .  AFWL  Eng  Model 

-  Cap  Model 

—  Elastic  Model 


Cap 


Elastic 


TABLE  5.  TIMING  COMPARISON  FOR  ONE-DIMENSIONAL  WAVE  PROPAGATION 


CPU  Seconds 


Calculation 

Model 

Total 

Plotting 

Grid  Setup 

Grid  Cycle 

ARA 

7114 

45 

2326 

4.22-  4.29 

Planar 

Cap 

1881 

45 

730 

0.45-  1.19 

(40  ms) 

AFWL  Eng. 

595 

45 

5 

0.48-  0.50 

Elastic 

135 

46 

4 

0.21 

ARA 

1247 

21 

5.19 

2.14-10.29 

Spherical 

Cap 

381 

24 

4.43 

0.42-  6.34 

(2  ms) 

AFWL  Eng 

151 

21 

5.13 

0.52-  0.54 

Elastic 

94 

21 

4.48 

0.25 

ARA 

5366 

52 

3.18 

3.38-15.60 

Cyl indrical 

Cap 

866 

54 

2.74 

0.11-  2.22 

(40  ms) 

AFWL  Eng 

438 

55 

2.96 

0.49 

Elastic 

146 

54 

2.77 

0.21 

totes:  (1) 

Ail  0  d  1  C  U  1  a 

ti  ons 

had  99  grid  points. 

98  zones 

(2) 

Computer  used  was 

;  a  Digital  VAX  11/750,  system  clo 

ck  resolution 

0.01  sec. 

(3) 

"Total"  CPU 

time 

includes 

the  entire 

cal cul ati on. 

(4) 

"Plotting" 

CPU  time  is  the  time  taken  to  generate 

the  time 

histories  and  grid  plots. 

and  is  not 

influenced  by 

model  type. 

(5)  "Grid  setup"  CPU  time  includes  dimensioning  zones,  numbering,  ct 
as  well  as  ini  tial  izatiori  of  zone  stress  to  account  tor  gravity. 


(6)  A  "grid  cycle"  is  the  time  reuuired  tor  one  calcu iatioria!  cycie 
through  all  99  grid  points. 


4.5  Two-Dimensional  Calculations 

4.5.1  Code  Description.  The  code  chosen  for  the  two-dimensional  wave 
propagation  calculations  was  STEALTH.  STEALTH  (Solids  and  Thermal  Hydraulics 
Codes  for  EPRI  Adapted  from  Lagrange  TOODY  and  HEMP)  is  a  general  purpose 
commercially  available  code  developed  by  Hofmann  (1978)  for  the  Electric  Power 
Research  Institute  (EPRI)  and  is  completely  documented  in  [Hofmann  (1981)].  The 
DOD  version  of  STEALTH  was  used  for  these  calculations  and  is  an  adaptation  of 
Version  4.1a.  The  principal  reasons  for  choosing  STEALTH  were: 

(1)  Prior  experience  with  the  code  for  blast  and  shock  problems 

(2)  Available  to  the  public  and  fully  documented 

(3)  Code  modularity  and  ease  of  modification. 

4.5.2  SEM-STEALTH/2D  Link.  One  of  the  common  sources  of  uncertainty  which 
arise  when  comparing  calculations  made  with  different  codes  but  ostensibly  the 
same  material  model  stems  from  small  coding  differences  between  the  model 
versions  actually  used.  This  uncertainty  can  be  eliminated  by  using  the  same 
physical  coding,  drawn  from  a  central  library,  for  each  code.  This  is  possible, 
and,  in  fact,  practical,  for  many  finite  difference  and  finite  element  codes 
because  the  required  material  model  formulation  is  often  code  Independent.  The 
inputs  are  a  strain  increment  tensor,  the  old  stress  tensor,  and  state  variables 
and  the  output  is  a  new  stress  tensor. 

Toward  this  goal,  the  Soil  Element  Model  (SEM)  may  be  viewed  as  a  library 
of  constitutive  models.  Linder  this  effort,  the  SEM  was  implemented  as  such  and 
linked  directly  with  STEALTH  2D  on  the  AFWL  CRAY.  The  way  in  which  this  was 
accomplished  is  outlined  in  Figure  38.  There  were  four  principal  places  where 
STEALTH  needed  to  be  modified  to  incorporate  the  SEM  models: 

(1)  The  main  program,  where  SEM  cornnon  blocks  holding  material  properties 
were  inserted  to  insure  contiguous  memory  locations.  The  actual 
link-up  between  the  programs  occurs  during  the  loader/linker  phase 
where  the  SEM  routines  are  made  available  to  STEALTH  as  a  binary 

1 Ibrary. 

(2)  The  material  input  phase,  where  the  Input  was  reconfigured  to  use  the 
SEM  parameter  Input  subroutine  (INPEOS),  In  this  way,  the  SEM  library 
of  model  parameters  could  be  read  by  STEALTH  directly.  Input  model 
parameters  are  also  echoed  using  the  SEM  routine  for  this  (OUTEOS). 
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Figure  36.  Soil  Element  Model  -  STEALTH  2D  Interface. 


(3)  The  problem  initialization  phase,  where  gravity  stresses  are  set  for 
each  zone.  Again,  a  SEM  routine  was  used  (GLOAD),  which  sets  gravity 
stresses  depending  on  the  appropriate  SEM  material  model. 

(4)  The  zone  processor  phase,  where  STEALTH  simply  calls  the  main  SEM 
model  subroutine  (MODEL)  which  in  turn  calls  the  appropriate  model. 
Some  of  the  SEM  boundary  load  application  routines  used  in  SNEAKY, 
e.g.,  Speicher-Brode  Nuclear  Overpressure  (SPBRODE),  may  also  be  used 
as  an  option. 

Thus,  all  the  models  implemented  in  the  SEM  are  now  implemented  in  STEALTH  2D. 
This  interface  was  used  for  the  ARA  model  two-dimensional  continuum  code 
checkout.  A  full  listing  of  the  updates  for  the  SEM-STEALTH  modification  Is 
provided  in  Appendix  B. 

4.5.3  Planar  (DIHEST)  Calculations.  Several  options  were  considered  while 
choosing  a  problem  for  the  two-dimensional  calculations,  including:  a  traveling 
nuclear  airblast  over  a  half-space  (HEST),  a  buried  explosive  sphere,  a  surface- 
flush  cylindrical  charge  (CIST),  and  a  vertically-oriented  planar  array  of 
burled  charges  (DIHEST).  All  would  have  been  adequate  choices  because  they  are 
all  used  by  the  Air  Force  in  simulating  one  aspect  or  another  of  nuclear  weapon 
ground  shock  effects,  and  all  produce  two-dimensional  wave  fields.  The  DIHEST 
was  chosen  for  the  test  case  because: 

(1)  It  produces  waves  which  immediately  interact  with  the  free  surface,  an 
important  aspect  of  many  two-dimensional  problems. 

(2)  It  has  a  relatively  simple  geometry. 

(3)  Many  DIHEST  experiments  have  been  performed  in  dry  alluvium  and  data 
is  available  for  comparison. 

DIHEST  events  are  commonly  used  to  simulate  upstream-induced  ground  shock 
effects.  They  have  also  been  used  to  simulate  earthquake-like  motions,  and  the 
geometry  of  the  problem  calculated  here  was  taken  from  the  SIMQUAKE  test  series 
[Higgins,  et  al.  (1983)].  Figure  39  shows  the  SIMQUAKE  II  experiment.  For 
these  calculations,  a  single  DIHEST  array  was  assumed.  The  calculational 
geometry  is  shown  in  Figure  40.  Note  that  the  grid  Is  inverted  with  the  ground 
surface  at  the  bottom  of  the  figure.  This  conforms  to  the  STEALTH  convention, 
and  is  the  most  consistent  way  of  visualizing  the  calculational  set-up. 
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Figure  39.  Simquake  II  Elevation  [Higgins,  et  al  (1983:2-4)]. 
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The  mesh  was  approximately  175  m  wide  by  133  m  deep,  and  had  2,925  grid 
points.  The  grid  was  large  enough  to  avoid  boundary  reflection  effects  (from 
the  top  and  right  side)  into  the  region  of  interest  within  the  time  of  interest. 

The  problem  was  run  to  0.5  second.  As  in  the  one-dimensional  calculations,  only 
quadratic  artificial  was  employed  (in  compression)  and  the  coefficient  of 
artificial  viscosity  (CQV)  was  set  equal  to  2.0. 

A  pressure  history  which  corresponded  to  a  SIMQUAKE-like  array  was  applied 
to  the  left  side  of  the  grid.  The  form  of  the  pressure  function  was  that  used 
by  Higgins,  Johnson  and  Triandafilidis  (1978): 

P(t)  =  P0(l-t/t0)e"a  t/l°  (0  <  t  <  tv)  (66) 

where  P  =  pressure  (Pa),  t  =  time  (sec),  P0  =  peak  pressure  (Pa),  t0  =  duration 
coefficient  (sec),  a  =  decay  coefficient,  and  tv  =  time  of  venting  (sec). 

The  peak  pressure,  P0,  is  a  function  of  explosive  type  and  was  taken  to  be 
15x10^  Pa,  which  was  used  In  SIMQUAKE.  The  duration  coefficient,  t0,  is  a 
function  of  explosive  loading  density  and  time  to  venting,  and  was  set  to  0.070 
sec.  The  decay  coefficient,  a,  is  also  a  function  of  explosive  loading  density 
as  well  as  soil  type  and  was  assumed  to  be  11.4  for  this  problem.  The  pressure 
history  and  associated  impulse  are  shown  in  Figure.  41.  ' 

The  target,  point  locations  were  chosen  to  provide  complete  coverage  of  a 
rectangular  area  within  100  m  laterally  and  60  m  vertically  of  trie  origin.  Many 
target  points  were  chosen  to  correspond  with  gage  locations  in  the  SIMQUAKE  I 

events.  Target  point  locations  are  indicated  in  Figure  40  and  listed  in  Table  6 
along  with  the  kind  of  plots  made  for  each.  A  view  of  the  distorted  grid  near 
the  source  is  shown  in  Figure  42. 

t 

4.5.4  Discussion  of  Two-Dimensional  Results.  The  ARA  model  was  used  to  I 

calculate  free-field  response  to  the  DIHEST  loading  out  to  500  ms.  The  same  ; 

model  parameters  used  in  the  one-dimensional  calculations  were  used  here,  except 
cohesion  was  set  to  zero  to  enhance  the  effect  of  the  free  surface.  , 

I 

Calculated  horizontal  and  vertical  velocities  are  shown  in  Figures  43  and 
44.  Several  points  may  be  made  upon  examination  of  these  figures:  ' 

(1)  From  the  times-of-arr i val  along  the  array  mid-depth  gage  line,  the  ! 

calculated  wavespeed  is  seen  to  be  321  m/s.  The  peak,  is  traveling  at  | 

230  m/s.  > 

(2)  The  peak  horizontal  velocity  along  this  sam-  gage  line  attenuates  as  j 

shown  in  Figure  45.  Near-surface  horizontal  velocity  attenuation  is 

also  shown.  ” 
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0 1  NEST  Horizontal  Velocity  Histories  Using  the  ARA  Model 


(3)  The  duration  of  the  calculation  does  not  permit  full  development  of  a 
negative  horizontal  velocity  phase. 

(4)  Horizontal  velocities  at  ranges  greater  than  50  m  tend  to  be  greater 
near  the  surface  than  at  array  mid-depth. 

(5)  The  spall  duration  can  be  observed  in  the  close-in,  near-surface 
vertical  velocity  traces  and  in  some  zones  near  the  edges  of  the 
source. 

In  order  to  better  understand  ARA  model  behavior  under  general  two- 
dimensional  stress  and  strain  paths,  calculated  octahedral  stress  paths  and 
pressure-volume  response  were  tracked  at  several  locations.  These  are  shown  in 
Figures  46  and  47,  respectively.  The  following  observations  can  be  made  from 
these  results: 

(1)  Loading  (even  directly  in  front  of  the  source)  rapidly  diverges  from 
uniaxial  strain  conditions,  and  the  stress  point  quickly  encounters 
tne  expansive  yield  surface. 

(2)  Due  to  geostatic  stresses,  loading  from  the  side  causes  an  Initial 
drop  in  shear  stress  (as  demonstrated  in  Section  4.2).  Very  near  the 
surface,  this  Initial  drop  is  not  observed  bpcause  loading  is  oriented 
more  toward  vertical. 

(3)  The  variability  in  stress  paths  is  due  to  the  many  different  kinds  of 
wa.es  emanating  from  the  source,  reflections  off  the  *ree  surface,  and 
numerical  oscillations  associated  with  the  grid. 

(4)  The  variability  in  volume  compressibility  response  is  due  to  the  same 
factors  mentioned  above  because  different  stress  paths  cause  variable 
activation  of  the  two  yield  surfaces  in  this  model.  Much  of  the 
differences  with  range  can  be  attributed  to  the  changing  peak  stress 
level . 

(5)  Spall  and  rejoin  behavior  is  very  important  near  the  surface  at 
close-in  ranges,  as  well  as  near  the  edges  of  the  source. 

The  array  mid-depth  target  point  at  61  m  range  is  a  good  example  of  the 
complex  loading  history  experienced  in  this  kind  of  event.  Figure  48  reviews 
the  stress  path  and  pressure-volume  response  for  this  location.  Corresponding 
points  on  the  two  curves  have  been  noted  and  the  complete  history  may  be 
followed  beginning  at  point  1.  Note  the  difference  in  effective  bulk  modulus 
depending  on  the  direction  of  the  stress  path  and  currently  activated  yield 
surface( s ) . 
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Figure  48.  Stress  and  Strain  History  at  Range  =  61m,  Depth  =  19m 


The  calculation  using  the  ARA  model  can  be  compared  to  SIMQUAKE  data, 
although  it  should  be  recalled  that  the  model  parameters  have  not  been  fit  with 
insitu  behavior  in  mind.  The  inability  of  the  lab-based  model  to  predict  insltu 
behavior  was  clearly  demonstrated  in  the  one-dimensional  calculations.  The 
behavior  there  was  generally  too  soft  (high  velocities)  and  too  slow  (late 
arrivals,  peaks).  Also  recall  that  SIMQUAKE  was  conducted  at  McCormick  Ranch, 
which  is  also  a  dry  alluvial  site,  but  McCormick  Ranch  sand  has  somewhat 
different  properties  than  does  CARES-OAV  alluvium. 

Figure  49  shows  comparisons  between  the  calculation  and  SIMQUAKE  II  data  at 
the  61  m  range.  The  overall  trends  appear  correct,  however,  the  timing  and 
magnitudes  are  not. 
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Figure  49.  Comparison  of  Velocity  Waveforms  Calculated  Using 

the  ARA  Model  with  Simquake  II  Data  at  the  61m  Range 
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5.0  CONCLUSIONS  AND  SUMMARY 


5. 1  ARA  Model  Assessment 

The  primary  goal  of  this  effort  has  been  achieved:  the  ARA  model  developed 
and  previously  tested  under  laboratory  stress  and  strain  paths  has  been  success¬ 
fully  Implemented  and  tested  in  dynamic  finite  difference  calculations.  Many 
aspects  of  the  model,  Including  numerical  errors,  cohesion,  work-softening,  and 
tensile  failure  Issues  have  been  adjusted  and/or  Improved. 

The  most  Important  aspect  of  the  ARA  three  invariant  model  is  that  it 
accounts  for  many  observed  features  of  soil  behavior  in  the  framework  of  a 
rigorous  plasticity-based  formulation.  In  this  respect.  It  Is  fundamentally 
different  from  most  models  currently  In  use  for  ground  shock  calculations. 
Intuitively,  one  would  expect  a  model  which  can  accurately  predict  the  labora¬ 
tory  response  of  soil  under  a  complete  suite  of  stress  and  strain  paths  will 
better  predict  Insltu  blast  loading  response.  The  calculations  performed  for 
this  study,  however,  do  not  entirely  support  intuition.  In  many  cases,  the  ARA 
model  did  not  produce  substantially  improved  wave  propagation  response.  Two 
questions  must  then  be  addressed:  (1)  Why  not?,  and  (2)  Can  this  observation 
be  extrapolated  to  calculations  of  ground  shock  In  general? 

One  of  the  aspects  of  soil  behavior  under  blast  loading  that  becomes 
Immediately  apparent  upon  studying  calculated  stress  paths  Is  that  a  given 
element  of  soil  spends  only  a  small  fraction  of  its  entire  stress  history  In 
initial  loading.  Typically,  the  soil  is  loaded,  unloads  in  pressure,  and 
quickly  approaches  yielding  in  shear.  Much  of  the  post  peak  response  is 
dominated  by  the  unload-reload  behavior  of  the  model.  Any  secondary  loading  is 
heavily  Influenced  by  this,  as  well.  In  many  cases,  the  soil  will  also  fail  in 
tension,  and  this  makes  post-spall  behavior  critical.  Thus,  a  fairly  general 
hypothesis  can  be  proposed:  any  new  model  which  concentrates  on  initial  loading 
and  plasticity  effects  related  to  Initial  loading,  but  reverts  to  simple  elastic 
behavior  upon  unloading  or  reloading  will  tend  to  produce  dynamic  waveforms  very 
similar  to  that  which  can  already  be  produced  with  simpler  models.  The  one  and 
two  dimensional  calculations  presented  here  provide  some  evidence  to  support 
this.  More  detailed  study  and  analysis,  Including  models  with  improved  unload- 
reload  behavior,  are  required  to  substantiate  the  hypothesis. 

89 


i 


K 


I 


t 


V. 


& 


k> 

k.« 


K 


s* 

>: 


i 


► 


Additionally,  the  treatment  of  tensile  failure,  where  the  material  In  a 
zone  can  become  distended  or  cracked.  Is  very  important  for  many  blast  loading 
geometries.  Directional  tensile  failure  is  particularly  tricky  to  Incorporate  in 
a  model  formulated  like  the  ARA  model,  because  it  is  not  clear  how  the  yield  and 
plastic  potentials  would  be  affected.  So,  much  of  the  similarity  between 
waveforms  produced  by  the  ARA  and  other  models  can  potentially  be  explained  by 
fundamentally  similar  treatment  of  unload-reload  and  tension  behavior. 

Given  the  similarity  of  results,  one  of  the  main  computational  deficiencies 
of  the  ARA  model  Is  Its  expense.  A  typical  wave  propagation  calculation  takes 
an  order  of  magnitude  longer  to  run  than  an  elastic  calculation,  seven  to  eight 
times  longer  than  an  AFWL  Engineering  model  calculation,  and  three  to  six  times 
longer  than  a  cap  model  calculation.  The  biggest  slowdown  occurs  when  strain 
subcycling  is  required  to  satisfy  the  consistency  condition  on  the  expansive 
yield  surface.  This  occurs  upon  Initial  departure  of  the  stress  point  from  the 
hydrostatic  axis,  which  usually  begins  during  geostatic  loading.  Although  the 
ARA  model  logic  Itself  Is  certainly  more  complicated  than  that  of  the  other 
models,  this  does  not  really  result  in  a  significant  time  penalty.  Improvements 
in  the  subcycling  technique  or  In  the  formulation  of  the  expansive  work  rela¬ 
tionship  (to  eliminate  the  singularity  at  zero  expansive  work)  would  be  most 
beneficial  to  overall  efficiency  of  the  model. 

Probably  the  most  important  aspect  of  the  ARA  model  calculations  is  the 
coupling  of  shear  and  volumetric  response  through  dilatation.  This  is  most 
observable  in  the  regions  of  highest  shear  such  as  those  adjacent  to  an  explo¬ 
sive  source.  However,  there  are  still  some  issues  which  relate  to  this  kind  of 
behavior  which  remain  to  be  resolved.  For  example,  the  ARA  model  has  a  separate 
plastic  potential  related  to  the  expansive  yield  surface  which  allows  for  direct 
control  of  dilatation.  But  how  is  this  plastic  potential  affected  by  unloading 
and  subsequent  reloading?  The  ARA  model  is  capable  of  dilation  while  unloading 
In  pressure  while  remaining  on  the  expansive  yield  surface.  To  determine 
whether  this  Is  physically  realistic  would  require  some  laboratory  testing.  The 
tendency  In  the  ARA  model  for  the  stress  point  to  move  very  rapidly  down  the 
failure  surface  toward  the  apex  seems  to  significantly  influence  calculated 
behavior.  In  fact,  this  Is  likely  the  biggest  reason  for  the  differences 
observed  between  the  ARA  and  cap  model  calculated  results.  This  phenomenon  and 
the  parameters  responsible  for  it  need  to  be  better  understood. 
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Many  questions  regarding  the  suitability  of  the  ARA  model  for  general  use 
In  ground  shock  and  soil-structure  interaction  problems  cannot  be  answered 
solely  on  the  basis  of  the  calculations  completed  under  this  effort.  Therefore, 
it  Is  appropriate  to  record  here  reconmendations  for  continued  Improvement  of 
the  ARA  model. 

5.2  Recommendations  for  Continued  Study 

The  following  is  a  prioritized  list  of  tasks  which  would  help  to  move  the 
ARA  model  into  the  mainstream  of  ground  shock  computing  and  improve  its  pre¬ 
dictive  capabilities: 

(1)  Improve  the  subcycling  procedures  for  correcting  yield  surface 
overshoot  due  to  both  extrapolation  of  afp/aWp  over  a  tlmestep  and 
elastic  punch-through.  The  goal  is  Improved  efficiency,  competitive 
with  currently  available  models  such  as  the  cap  model. 

(2)  Re-examine  the  work-softening  aspect  of  the  ARA  model  to  determine  if 
it  can  be  successfully  incorporated  in  two-dimensional  calculations 
without  violating  uniqueness. 

(3)  Reformulate  the  compressive  yield  function  (f^-)  into  a  dimensionless 
form.  This  will  eliminate  the  very  large  numbers  encountered  with 
some  sets  of  stress  units  (such  as  Pascals)  and  allow  the  model  to  be 
used  on  smaller  computers  in  these  units. 

(4)  Fit  the  ARA  model  parameters  to  stress-strain  data  which  is  most 
representative  of  Insitu  behavior  for  the  specific  site  of  each  field 
test  and  rerun  the  dynamic  calculations.  Evaluate  the  ARA  model  in 
Its  ability  to  predict  blast-loading  events. 

(5)  Formulate  a  rezonlng  strategy  which  accounts  for  a  redistribution  of 
the  state  parameters  In  the  ARA  model.  Rezonlng  Is  a  nasty  fact  of 
life  and  Is  used  often  In  severe-envlronment  calculations.  It 
Involves  the  smearing  and  lumping-together  of  greatly  compressed  or 
distorted  zones  to  allow  the  calculation  to  proceed  at  a  reasonable 
tlmestep.  The  several  state  variables  used  by  the  ARA  model  must  be 
consistently  transported  from  old  zones  to  new  zones.  Note  that  this 
is  not  a  problem  unique  to  the  ARA  model  and  must  be  dealt  with  when 
using  any  constitutive  model  with  history-dependent  parameters. 

(6)  Implement  the  unloading-reloading  hysteresis  coding  which  has  been 
formulated  based  on  a  Masing-like  relationship  for  bulk  moduli.  Test 
out  the  formulation  for  stability  in  a  finite  difference  calculation. 
If  this  approach  fails  to  produce  improved  unload-reloading  cyclic 
behavior,  consider  Incorporating  kinematic  hardening  relationships 
Into  the  model. 

(7)  Perform  sensitivity  studies  to  evaluate  the  effect  of  various  model 
parameters  on  calculated  waveform  features.  This  is  a  critical  and 
necessary  step  In  the  acceptance  of  the  ARA  model  as  a  viable  ground 
shock  model . 


(8)  Couple  the  ARA  model  with  a  good  tensile  failure  model  which  allows 
directional  cracking.  This  appears  to  be  necessary  for  many  kinds  of 
loading.  For  example,  In  an  axlsymetrlc  explosion  geometry,  It  is 
desirable  to  allow  the  material  to  separate  In  the  hoop  direction 
while  still  transmitting  stress  in  the  radial  direction. 

(9)  More  fully  test  out  the  high  pressure-temperature  equation  of  state 
developed  In  this  effort  In  appropriate  calculations,  perhaps 
culminating  in  a  two-dimensional  nuclear  source  energy-deposition 
calculation. 

(10)  Upon  completion  (or  at  least  resolution)  of  Items  1-9,  use  the  ARA 
model  to  calculate  a  more  complicated  blast  loading  problem  of  current 
interest  to  the  Air  Force.  Two  suggestions  at  this  point  In  time 
would  be  a  combined  HEST-DIHEST  simulation  or  an  NSS  event.  However, 
there  may  be  more  appropriate  choices  when  the  time  comes  to  seriously 
apply  the  ARA  model . 

(11)  Because  the  ARA  model  can  accurately  model  shear-volume  response 
coupling.  It  Is  a  good  candidate  for  the  grain  matrix  in  a  coupled 
fluid-mechanical  calculation.  Incorporating  this  model  Into  a  dynamic 
fluid-mechanical  code  should  be  considered. 
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APPENDIX  A 

LISTING  OF  THE  ARA  MODEL 


applied  research  associates  three-invariant  soil  model 
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=  OR  COMPLEX  DYNAMIC  LOADINGS:  Ce*elD;^er,:  of  a  Tnree 
I^va'iant  Constitutive  Mod*  1  "  .  A  F  OSR-T  R- 8  5-XXXX  ,  Aac'led 
Research  Assoc -fates  ,  Inc  .  .  Albuque-qu*.  NM.  Aprfi  1S35. 
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COMMON  /  A  _  ^  MOD/  I  ECS  .  ;  UN  I  T3 atmc,  CCNV  ,  RhCREF  .  “T  VPE  ,  GR  .  GRP  , 
85P(3).VSTNl3).8MODN.  uMOON  . STA  T£ .  ‘  C  ) 

COMMON  /CCNSTAn  */wGNthC. TCC ! H U . S(JR  T  2 , SORT  3 , SQRT& , SQRT  2 3 , 

TwOP! .eiGPCS.BIGNEG 
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CLV,  RDS,  XACGR,  NMP,  NMAT  ,  ITPS.  IFF  T.  PI  .  I  £  ,  *XE.  ~EG 
AV’Py3) .SMLNUM.BIGNJM.MOOPRT 

DIMENSION  S:GC(6),S!GI(e;,DEPO(6).OEPI(6),S!G(6).CEP(6) .SDEV ;S) . 
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81(6), 82(6). ACE  ( 6,6).ACEP(6,6),DSIG(6). DFPPDS(6 )  . 


/POOP ' 
/SN  E  A  X 


.N  .N  - 


*  OG P OS ( 6 ) ,F8AK(6 . 2 ) . DS0(2. 2 ) , A0AR(2 ,6 ) , ALE AP( 2. 6)  , 

*  Ai_6AR](2,S),B6aR(2,6),GBAR(6.2),  DLAM( 2  )  .  AC  EG  (6  }  . 

*  TCEP(  6 , 2  )  .  OFEP(2) . AMCOLI6)  , STATSAVC  8) 

REAL  J2.J3 

OENCONV  AND  ENGCCNV  ARE  CONVERSION  FACTORS  TO  CONVERT  DENSITY  TO  GM/CM 
ANO  ENERGY  DENSITY; ENERGV  PER  -NIT  MASS)  TO  TERG/GM  IN  THE  HPEOS. 

DIMENSION  OENCONV ( 8 ) . EN3CCNV ( 9 ) 

data;  GENCCN  V(  r),I-',S)/1  .  G  E  -  3  .  '  .  C  ,0.0.0.0.00.’.  E-S,0.0.  ’  .  0/ 
DATA(ENGCONN,(  I  ;  .  1  =  1  ,  8  )  /  1  .  C  E  8  . GE1 2  . 0.0.0.  0,2.0.  ’  .  0E8  .  0 . 0 ,  1  .0/ 


DMN  IS  A  MINIMUM  DENSITY  (G/CC)  FOR  TH  E  HPEOS 

da ’a  /c . : •/ 


:’S’=MC>(T  c- A,r*  IONS 


■ . ;o:~  ■ =so;t» af ex 

=  D  r  *  3D  ;  ~  ^  too  t  1=3  -ar-C(SCCt  ;  ■  *  D  *  3  •-CC-**! 

= 2C'  (  AhD  ;  =  a  •  V.;  -  *2  *  •  a  wC  /  I  ACC  »  A  *v.y  1  *  * \-L2-) 
'I-  S'-SEE  =av«  3-0333  atm.  d  _■  > 


-  ’  -  :«t  i  hhj  ;■  3  ,  e  =  A C  * 

if  •  -1E*A^P/A-VD  }-£/•-•.  -  A  0  •  A  -  /  A  ’ 

I  ,  AhP  ■.  ArvP  i  A  VA  ,  •  a*D  ,  1  £  -  •  5,  ) 


AA/A  -M  D»  (  AO*  ExP  I  -  a  c  *  A*.c/A"MD  )  -  Ac  •  EaP(  -  AS  *  awP/A’MD  )  ) 

F  •  ,  jD  .  ^3  ,  SGC  *  )  =  ».CNth  D*  ,  S0R*2  3*  -  2**  3  .  E-2*  a£v*--  S0R’2*. 

,  2. SO- ’,**2 

c  2 : J  2 . . 3 . SCO  *  )  *  (  ■  /SORT c  5 . * J  2  )  -  3*  t E v* -  3/ I s;p*  2* J  2  *  *2 , ) 

*  ;  ’  /  C  ;  SCO."  )  -amv  /A’MC  ) 

p  3  v  2  .  SOC  *  '  =  3  *  A  E  v/  I  SCR  *  2  * .'  0  )  *  I  ' ,  0  I  SOD  *  )  •»*»  ,/jtmd  •; 

-Y  ,  A  E  ’  i  I  ,  '  ■■  '  =  A  E  *  A  2  ;  •  A  r  M  c  '  *A**-*D,3_I-  .A’V  *  *  , 

o2;j2.-3;=,'.o;p'';E.*.2;-'3*ae-*.:/\SOF"-*-2**2,  'a-md 
32(-,2,*J*a£y/i  3v;.*2*_-2*atmo> 

SIGN  .-MAN  G  F 

CA-.L  MSCA-R(SIGC  .  6  .  ’  . -1  .  2  ,S!G!  ) 
tall  mscalrcs is: ; a  > , 3 . • . sqr~2. sigi (a ) ) 

CALL  msc AL  RfOERC , 5.  1 , -  1  . 3 . DEP!  ) 

CALL  MSCALR(S:G! ,fc, 1 . 1 .C.SIG; 

INITIALIZATIONS 


mt y°E-0 
OWC-O . 0 
DNP=D . 0 
I  SUB-0 
NSUB-0 


5v' 


5 


5 

L  .V 


i 


r: 


/ 

V 


V 

V. 


0115 

OVSMOOT  *0 . 0 

0116 

10  DO  11  1*1 ,NSVAR( 10) 

0  117 

11  S"ATSAV( I)=STATE( I  ) 

0  118 

c 

0  l  '9 

c 

RESTORE  STATE  VARIABLES 

0120 

c 

0  121 

20  ASOCT*STATE( 1 ) 

C  '  2  2 

A«C*S’A-rE(2) 

2  123 

AwP=STATE( 3) 

3'  24 

A . 2  =ST  AT £ ( 4 ) 

0125 

ASPALL*STATE(5 ) 

2  126 

ACONFP*STATE( 6 ) 

2  1  2T 

A=V=STATE(7) 

2  1  2? 

A3=ENG=StATE(83 

:  ’  29 

a:vap  =  s''AtE(S; 

2  3  C 

c 

•  •  •}  ' 

?AvE  INI'IA.  ;~f£33  o;  ; 

.  :  3 

: :  5 ; : -sav=scc* 

*  -1 

*""3AV  =  ~c :  ” 

'  i  5 

-  2  -  3  A  v  *  F  2  0  ;  "C'C  T  ,  ”CC  T  ; 

7  J  5 

cC2PSAV*=C2P(AL.C; 

:  •  ?  7 

cdpSAV*-pp; SOCT , TCCT  .  CC 

;  "  5 

-  4 2 P 3 A V =  =  P2 P; A«*P; 

'  *  39 

' 

.  1  4  C 

z 

<*:  d;fv  incoming  stress 

.  4  ‘ 

' 

;  ■  ;  3wE~rac~  out  vapor 

:  *  5: 

l 

•  3 

IE  CO  37  1=1.5 

1  4  ‘ 


149 


0  15  1 
C  ‘  5  2 
2  ’  5  3 


0 


0‘55 
0  156 
0  ’  57 
2  •  58 
0 '59 
0  ■  50 
e 1 

2  '52 

:  •  53 

0 ' 64 
0  165 
0  ‘66 
0167 
C  i6e 
0169 
0170 
0171 


.1  \  3  <•  ’  r 


5 .  3 1 1 1 =5 :  5 : ;  ■  )-apvap*a*<;:ol  c  r ) 
24  LL  MSCALRI 3131 . 6 , 1 . '  2  ,  S  '  G  ) 


initialize  energy  related  variables 


ET  A= 1  ./( 1  • -AEV) 
D;N3«RHOREc*ETA 
ENGN*Z I E/RRGREF 
DENG=ENGN-ASPENG 
ASPENG=ENGN 
kE*ENGN/aE^El I 
E  F  AC  = 1  . 0-RE 

;p(A l-E.NE. I . 0j£RAC*l . 0 
OELMU  EN “ACTE* DENG* AH E 
Ic<  EFAC . LT . 0 . 0 )DELMLEN=0 . 0 


CALCULATE  MECHANICAL(DEVM)  AND  ENERGY(DEVE)  VOLUME  strain  INCREMENTS 


0Evm«0EPI( 1 )+0EPI(2)*DEPI(3) 
DEVE«DELMUEN/ETA**2 


INCREMENT  TOTAL  VOLUME  STRAIN 


OEV-DEVM+DEVE 

AEV=AEV+DEV 


C 

C 

C 


ADD  ENERGY-RELATED  STRAIN  INCREMENT  TO  ORIGINAL  STRAIN  INCREMENT  TENSOR 


98 


I 

I 


1 

i 

I 

I 


0  172 
0  1  73 
0  1  74 
0  1  7S 
0  1  76 
3  177 
3’  78 
0179 
0180 
3  •  8  1 
3  182 
3  '  e3 
C  '84 
2'85 
3  ’  66 
C  '  6* 
3  93 
3  '  93 


i  s  3 
:  ■  3« 
:  ■  35 


■  3  9 


DEPI  ( 1 )-DE»I ( ’  )*DEVE»WONTHO 
0EPI(2)»D£PI  (2)OEVE*WONTHO 
DEPi ( 3) -OEP I ( 3)-*-D£VE*WONTHD 
C 

C  SET  INITIAL  STRAIN  INCREMENT 

C 

CALL  MSCALfi(OEPI .6. 1 . 1 .O.OEP) 

C 

000’  IP( ISTOP.EQ. 1 )G0  TO  9300 
C  C-EC1  FOR  SPALLED  CONDITION 

C  ASPAlL  TRACKS  ’HE  EXPANSION  VOLUME  S’RAIN  WHEN  T-E  MA^ER I AL  EXCEEDS  ThE 

C  SPALL  ^IVIT  (APEX)  an 0  BECOMES  £XT£NOED.  WHEN  ASPALL  IS  ZERO  (OR 

C  POSITIVE).  MATERIAL  IS  COMPRESSING.  WHEN  ASPALL  IS  NEGATIVE. 

C  IS  SPAlwED.  ON  REjOIN.  wSE  CNlV  PORTION  OF  STRAIN  increment. 


•  z  '  ^  c  '  \.'.n  * r>  c  - 

4  A ; P; AV  =AGPA . _ 

4  3  i  c  ;  ASP  4  L  i. .  GE  .  3  3 9; 

C'E  V  I  =  1.  OEP I  (  ’  )  *CE  p  I  ;  3  I  *  0  E  R I  (  3  ;  ) 
ASPA^^-ASPA^L-’-  jEv  I 
IFvA3PAl>  l£  0  3)33  "3‘  6CC3 
Epa*to*aspalL/DEVI 
DO  45  I = • . 3 

44  3EP!'.  I)  =  E04*:c*35DI(:) 

=  5  3 £p; I ;=DEP; ( I  ) 

43  =  A'.L  =  3  .  3 


4  _■  2  C 


34 

0  5 
3  5 


3  2  C  8 
3  23  9 


32  '  3 
3  2  ’  4 
r  2  1  5 
32'6 
32  *  7 
32:6 
32  19 
3  2  2  0 
C  2  2  1 
0222 
3223 
0224 
0225 
0226 
3227 
0228 


S*R4IN  increment  «ITH  FP2p=0.0  (Al 2*0.0) 

'3  3  3  4  ..  ALA  IN S IS. SOEV . 3C3C . J 2 , 3 3 , S3C* , ~CCT . COS3W; 

C4_  ARASTA"-(S:G.  ACE) 

1 ’ 0  4PCP=FCP(SCCT , TOCT) 
a=C2P«c;2A,  4WC  ) 

!CC (4CC2P-AFCP) .GT. TCL2JGC  TO  303 
C 

2"  3  -3  4  MADC‘(S  IG,  5  ,  '  .  AMCCL  .  2  .  C  .  (  2  .  *(  4R-1  .  )*SCC *  >  ,  DRCOOS  ) 
C<_.  «;ci  R  ;  3 c  -3  p  3  S  .  5  .  ’  ,  '  .  3  .  30333) 

D  =  C2PDw*A  rMC/  (ALL  *4  PC  )  *  ;  4  CC  2  P/A*m;,  •  «  2  )  **  (  1  .  -  4  PC) 
HC=2*AFCP 

3'S0  (  1  .  1  )  *  0FC2PDW*HC 

C  A w  L  MTMULT( OFCPDS. 6 ,  1 , ACE , 6 . 6 ,  1  . 3 , A  1  ) 

CALL  M«LLT( A  1  .  ! , 6 , DGCDS, 6,1,1.0,RL’’) 

Ri_  1  1  *Rl  1  1  »DSQ(  1  ,  1  ) 

RL 1  *  1 .0/RL1 1 

CALL  MSCALR( A1 . 1 . 6 . RL 1 .8 1 ) 

CALL  mm'jlt ( 9  1  ,  ’  ,  6  ,  OEP  .  6  .  1  .  1  .  0  ,  DLAMC) 

DWC*HC»DLAMC 

IF(DWC.GT.O.O)  GO  tC  430 

c 

300  C4LL  MMULT( ACE . 6 . 6 . OEP. 6  .  1  ,  1  . 0 . DSIG) 

MTYPE-0 
GO  TO  500 
C 


400  AHC*AWC»OWC 


0229 
0230 
023  1 
0232 
0233 
0234 
C  2  3  5 
0236 

02  33 
22  39 
0240 
024  t 
024  2 
2243 
0  244 

-  -  *  c 

0  2  4  s! 


:  -  4  o 

2  2  4  9 

2  2  5' 
2  2  2 

2  2  54 
2  3  5  5 
22  5  5 


2  2  5  5 
2  2  6  0 
2  2  6  " 
0262 
0253 
J  2  6  4 

nej 

2255 
2  2  5  7 
02  59 
0269 

j  2  w 
2  2  7  * 
2  2  "  2 
2273 
0274 
0275 
0276 
0277 
0278 
0279 
0  2  8  0 
0281 
0282 
0283 
0284 
0285 


CALL  MADD( ARAEC , 6 . 1 , 0GC3S.  1.0.  DL AMC . ARAEC ) 

CALL  *i«LILT(ACE.6.6.DSCDS.6  .  '  .  1  .0  ,  AC  EG) 

CALL  **M'JLT  (  AC  EG ,  6  .I.B’.’.S.I.O.ACEP) 

CALL  MA20( ACE . 6 , 6 , ACEP .  1 .0.-1  .  0 . *CE P ) 

CALL  NMULT( ACE®, 6 , 6 . D£ P , 6 . 1 . 1 . 0 , DS I G ) 

“TYPE- 1 

500  C  A  _  L  **ADO(S!G,  6.  1  ,  OS  I G ,  1  .  .  1  .  ,S!G) 

CALL  ARAINV;SIG  ,  SOEV.  SCOP  .  J2  ,  J3  .  SOCT  ,  TOC7- .  COS3W) 
:cvJ2.LE  A  T “0/ 1  .  E  '  0  )  GO  TC  2000 

F I RST  ACTIVATION  CP  E*PANSIVE  *OCE 


522  A. 


:P7Lirrp(3CCT,''OCT.C053*»I 
'  4  '  4  -  E  p  =  2  r  ■  /  .  4  4  •  ;  A  0  -  4  - 


!  :  '  '59  .  7 ”  '22I3”2:556 

:  /  a  3*  a.cgc  ex? ' -Ae»7R : 4_ \-s: 

:  P  435(  A*.E«*-7«  I  A.)  .  S“_NL«  ■  ;GC 

-  R  ;  J.  *  AN  £W 
30  TC  555 
: * :  -  anew*  4t y  ;■ 

A-.  -  -4i*P*J)A.P 


' j  l 

570 


3£'.p®AL  CAS; 


.4-2=1  Cl 


755 


7  ■  C 


720 


800 


900 


2  4, _ A4AINV,3!G.SD£v.Ei-2r,j2,J2,3CCT,''C!C''.  0  C  S  7  w ) 

ca_ .  akastat(s;s. ace ) 

IP  MA'ERIAL  HAS  FUllY  *ELTED.  USE  ELASTIC  *CDULII 
:  P,  ECAC.  (.  =  .  2  .  C  )GO  TO  1300 

COMPARE  CURREN”  STRESS  POINT  POSITION  WITH  LOCATION  OP  EXPANSIVE  AN 0 
C cvppess ! VE  SURFACES  *C  2EtSRm:nE  '23SIS. E  rOOES  0r  YIE.CIVS 

AFC?=PCP( SCOT,* OCT) 

ACC2P«PC2D( AWC) 

A4;o=poo^3oo*,70CT,CC33w) 

APD20,FP20( awp ) 

IP'.  AWP  .  GT  .  AWPPK)AFP2P*4PP2P*AWSOPT*  =  P2P{  AWPPK)*(  1  .  0-AW3OPT) 

AE*A2  «  AET  A  2  0 ♦ ASG*  AP  P  2  ? 

I  P C ( AFC2P-AFCP) . GT . TQL2 )  GO  TO  900 

I F( ( AFP2P-AFPP) .GT . TOL2 )  GO  TO  1100 
GO  T0  1000 

1F( ( AFP2P-AFPP) .GT . TOL2)  GO  TO  1300 
GO  TO  1200 

STRESS  POINT  IS  AT  THE  CORNER  OF  THE  COMPRESSIVE  AND  EXPANSIVE  SURFACES 


y.f, 


1 


I 


3 


v 


El  * 

> 


f  v 


ft 


H 


V 

V 
> 


i 


1 2  6 


l  -  - 

o  i:e 
3  3  29 
03  30 
3  3  3' 
0  332 
0333 
033* 
0  335 
0336 
0  3  37 
0338 
0339 
03*0 
034  1 
03*2 


0286 

1  coo 

CALL  MADD(SIG. 6 . 1 . AMCCL ,2.C.(2.*LAP-,.)*S0CT),0FCPC3 

0287 

CALL  MSCALR( OFCPDS .6.1.10. OGCOS ) 

0288 

DFC2PO»*A'rMO/(ACC*APC)»( ACC2P/*T«C*»2)»*( i . -APC) 

0289 

“C “2*AFCP 

0290 

RF1-F1(J2. J3.SOCT) 

0091 

RF;-e-2(jj,  J3.SOCT) 

C  2  9  2 

RF3«F3(J2,SOCT) 

0293 

RG 1 *G 1 ( AST  A  2 , SOCT ) 

C  2  9* 

RG2  «G2 ( J  2 , J  3  ) 

0295 

R33=G3( J2) 

0296 

CALL  ARAFUNC ( 2  2. RF  1 . RF2, RF3. SOEV, SCOF,  DFPPOS ) 

0297 

CALL  ARAFUNC ( J2 . RG 1  , RG2 .  RG3  .  S  OEv  , SC OF . DGPDS  ) 

0298 

DFP2P«DCP2PD*'  AWD) 

0293 

:r(Awp.GT.  dc  p2  p  =*  op  0  r  d  •  a 

:  :  o 

- p  =  - RG  '  *  3  *SrC*»2*FG2*-2-3* “o  3* . 3 

“in* 

00  1331  I-'.S 

■  w  "  *1 

y : , '  ■  -  =•:  3 ;  ? ; 

:03 

7 : i ;  ( ;  ■-  --53  =  ,!; 

;  *?:a 

ssar- : . ■  =:c:2S( 

*  iq5 

ip- 

g?a  =  '  :-p?s ;: ) 

:  ■ :  e 

OSC'v  ,*0 

-  ?  -  7 

03 0( 1 ,2;-0.0 

nee 

333(2.  -)*3.2 

:  <09 

- -CL  2 . 2 ; =DF  ?2P*-P 

■  3  '  Z 

: a  _ w  v* V  -  z 7  ^  t  ^  z ,  ^  .  b  ,  ■  .  0 .  a  =»  a k  \ 

’  j  * 

0  *-L  VI*..  ,  i  =  AR,  2 . 6  ,  GBAR  ,  6 , 2  .  i  .  O.A-BAR) 

2 

CALL  MA„j;A.tAfc',0.2.OSv.  *  .  .  *  .  ,  A  L  6  A  P ; 

:  j  ■  3 

*  i  •  i 

a.-  4??i  ,s  a_  ’  ’  ,  C  '  *  #  .  5  a  >  '  f 

‘  •  5 

U  •  6 

c 

v  ;  0.  £  rr*£9v.'.i';.\  .CORNER  L\_i; - ao.je  yc'-;; 

2  •  ? 

0  1  ^ 

:e{  ALBAR.;  -  ,  •  ;  -E.3  CR  *-SAR;2,2).u£.C.  .CR.A_;aP  - 

:  ;  ’  fi 

*  •“  ~  ?  ‘39 

C  :  '9 

'C  ■  ■ 

3  ON T INI E 

0323 

C////Z 

S'CP  CALCULA'ICN  -  NGNUMCuE  MODE  SELECTION  DETECTED 

:  32  1 

ISTOP-1 

3  2  22 

GO  TO  0301 

:  323 

'  3  •  9 

CALL  “«ULt( A8AR, 2 . 6 , PEP. 6 .  1.  1  . 0 . 0CEP  1 

:  324 

'323 

CALL  POLAR'  , A .3  A  R (  1  .  1  1.A_SaR(2,  1  ,  . AL '  . RiO) 

'  . A_5*“(  '  . 

1 .  :cep(  • ) . 
.  RPQ  .  ANC. 


2)  .  A.S4R; 2 . 2 } . Al? . RSI ) 
DEEP' 2  > . OF  £ . 0MEG ) 

OMEG . LT . PS  I ,GO  TO  ’030 


. -(rW0PI/4.)  .ANC.  0**£3. 


3¥EG  .  .  v  TWO? 


E  .  R-0,'30 
/’2  )}GC 


:«_l  :o_a  = 

CALL,  POLAR 
1  C  29  !  F ( OMEG . GT 
.’  c(  OMEG . GT 
I  f  i  0*E  G  .  GE 
GO  TC  1  3  0  C 

3/////ECP  MODE  ACTIVE 
1  0  3C  DLAMC»0PE**L2/ALCAP*SIN( PSI-OMES) 
CLAMP»DFE*AL  */ALCAP»SIN( OMEG-RHO) 
OWC»Ol*MC*hC 

OwP«OLA>*P*mP 

ALFAC* 1 . C/ALCAP 

AL8  AR I ( 1  ,  1  >ALe*P(2.2) 

ALSAPI ( 1 , 2)«-ALBAR( 1,2) 

A -BARI (2 .  1  )--ALBAR(2.  1 ) 

ALB*RI  (2 . 2)«ALBAR( 1,1) 

CALL  MSCALR( AlBARI .2.2.ALFAC.ALBARI) 

CALL  MMULT( ALBARI . 2, 2. A8AR. 2 ,6,1.0,  BB AR ) 


:*o 


rO  1 0  5 1 


101 


C  34  3 
3  34ft 
034  5 
C  34  6 
0347 
0  348 
0  349 
0  350 
035  ' 
0  352 
0  3  5  3 
o  ?54 
0355 
-  ->  e  c 

0  3  5  ’ 

■  ;  c  a 
■  5  i 


-  ■  j 

.  3  74 

■  37c 
'  '  ’5 
'377 

3  3’8 

0  3  7  9 
0  380 


~  .  e  ? 

0  38  4 
3  3  95 
; :  a5 
3  3  9  7 

0  3  88 
0  389 
C  390 
039  ’ 
0392 
0393 
0394 
C  395 
0396 
0397 


GO  *  0  "800 
C/////6C  •'ODE  ACTIVE 
1040  DLA*C=0*£P( l )/4LB4R( 1 , 1 ) 

C'wCpQlAMC'HC 

cal'.  •*7»10LT(OFCPOS,  6  .1,  ACE, 6, 6,1  0  .  A  1  ) 
A'.pAC-  ’  .  3/A. 9  AR(  1  ,  ’  ) 

CALL  MSCALR( A1  .  1 , 6 , ALFAC. 8 1  ) 

30  T0  1600 

C/////EP  MODE  ACTIVE 
'.'50  DlA»p»DF£P(2)/AlBAR(2.2; 

OwP*DlAMP»HP 

CA.l  *«tMULT(  OFPPOS  ,6,  1  , ACE, 6, 6,  1  0,A2) 
AL  F  AC  *  1  0/Al8AR(2.2,' 

Ca__  *,0;ALp,i2.',6.4.ci0.°2) 


*  -  -  '•ao'C^s:  3.6,  • .  a*ccl.  '•  o.  ;2.  ".as- 
o-scalp  :ccoos , 6  .  ’ .  • .  o , occ: 3 ; 

:  =  ':  =  ;«--a-v;.  ,acc*apc)4-.a=c:-/a'''0**2)**(  -a 

'2  *Ar‘A 

-ic »ofc : po*.*p|c 

-A..  M-MJLT • 7FCPOS. 6 .1.ACE.6.6..0.A1) 

.  A  .  .  '  ,;,3'j.'C3,‘,  ’  .  '  1  .  a  _  i  •  \ 

:  _  ■  ■•030'',''' 

:.-s;=  :,-c_  • 

a..  v-.a.a,a.',5,a.'a;,3 
4  -  -  M v  t  ,  0’3  '  0  .  . .  ap*C  ’ 


o.  *o 

S  "  R  ESS  P  C  :  N  *  CN  EXPANSIVE  v;cLc  3l>Rp  AC  E 
•200  RF1«F1(J2.„3.SCC*) 

pc I*p; ; . 2 , j  3 , SOC” ) 

Rp3«=3' J2.SC0T) 

SO-  * C-  ,  A E  T  A 0 , 3 0 C ’’  ) 

»':-2*32  t  .  2  .  J  2  I 
R  3  *  3  3  :  .  2  ) 

CALL  ARAFUNC(,J2.»F1  ,RF2.RF3,30EV.SCOP,OPPPCS) 
CA.l  APAF  JNC  C  J  2  ,  RSI  .  RG2  .  RG3 . 3  0EV  .  5CC-C  ,  OGPDS  ) 
OF  =  2P«C  =  P2PD*. AwP) 

IP;  AWP.GT  .  AWPPX/'0CP2P»0FP2P*ARSCFT 
hP--RG1*3 . »SOCT  *  2*RG2  * J  2  -3*RG3*J3 
DSQ( 2 . 2)-OFP2P»HP 

CALL  MTWULT  (  DPPPDS.  6  .  1  ,  ACE  .  6 . 5  .  1  .  0  ,  A2  ) 

CALL  **MULT(  A2  .1.8.  OGPDS,  6.1.1.0.AL22) 

AL22«AL2  2  *OSQ( 2,2) 

ALFAC-1 . 0/AL22 

CALL  MSCALRt A2 , ’ , 6 . ALCAC . B2 ) 

CALL  MMULT (82, 1 . 6 . OEP , 6 . 1 . 1 . 0 , OLAMP) 

OWP«OLAMP»HP 

I F ( DWP . GT . 0 . 0 )  GO  TO  1700 


C 

C 


Elastic  stress  increment 


L*.t 


i 


l 


3 


i 


i 


f  * 

s 


K 


i 


a 


f. 


y 

y 


04  C3 

c 

040  ! 

'  300 

CALL  MMULT( ACE. 6 , 6 , 0EP.6 , 1 , 1 . . OSIG) 

0402 

1  304 

CALL  MADDCSIG.6. 1  . DS IG . 1  .  .  1 . , S 1G ) 

0403 

MType-0 

0404 

’305 

CALL  ARAINV(SIG.SD£V.SCOF. J2, J 3 , SOCT , TOCT . COS3W ) 

0405 

C 

0406 

c 

IP  ELAST!0  TRIAL  EXCEEDS  THE  EXPANSIVE  v;elD  SURFAC 

04  0  7 

c 

APEX.  CORRECT  STRESS  POINT  BACK  AT  CONSTANT  SCCT 

0408 

p 

0  409 

:  305 

IF( DWP . LT .0.0  . CR .  owe . LT . o . C )GO  TO  1319 

04  »  0 

I c ( £F AC  .lE.O.O )GO  TO  ’319 

04  1  1 

N*«0 

04  1  2 

EFCP«FCP(SOCT .TOCT) 

04  •  3 

£CC  2  P«CC2P ; AWC ) 

-  *•  **• 

:cpd  =  coD(SC01’,'‘OC*,0:OS2»'3 

34  1  S 

£C£>2P  =  PP2b;  AWP  ; 

•  e 

r  t  p  ^  15  =  E  r  ^  C  r  *  4  *.3  ^  r  ~  '  p  2  -  ,  4  "  c  k  •  •  '  ♦  - 

‘  X  ~ 

•  *  3 

•.  -  =  * 

*  x  '  * 

‘3-3 

: - t  V EFP?P-SPPP> . 3’ . TO_2 }GO  "O  1309 

-  4  J  C 

I  F  ,  N  «  .  E  0 . 0  :  N  K  =  0 

O'  4  *>  ' 

:s;nk.ec  ’iNK=3 

C  i  t  t 

3  0? 

O'ONT  ;nue 

:423 

: =  c  n  k  .  e  o . :  i  c-  o  ■ :  319 

:x2c 

:c((({3::-!i'*5!3(0)^:3(3):/3.o;.j-?«;  4*«c/*c  3 

,*  £  J  c 

3  ' 

-i':?=*o:'p'.sc:*.App2P)/'o:* 

e  4 :  * 

:  0  '  2  •  j  :  -  ■ .  3 

34:^ 

":-tv( :  ;»s oev;  : . * B a t  1 0 

*  4  r.  m 

3  ’  * 

•:  4 : ; 

’  *  ‘  - 

ar  a  : '  v ,  3  :g  .  so-Ev .  soo r . 0 .  v  ?  .  sccT  .  Tc  rT .  o:  3  3  a  , 

C  4  30 

*  0  :  3  ‘  s 

0  4  3’ 

0 

, 

C 4  30 

V 

IP  STRESS  POINT  EXCEEDS  EIT-ER  v:E-D  SLRPACE  BEXDNO 

04  ?  3 

c 

SUBCYCLE  5  a  5E  0  ON  The  pa t;o  F‘/f‘.  PRE-EXISTING  O'. 

04  34 

c 

CANNOT  6 E  CORRECTED  AND  IS  EL6',‘EaCtE0  OL'  . 

04  35 

' 

04  36 

'  3 ' 5 

2  p ( N3UB  ED .  C)GO  TO  3  316 

043*7 

GD  ~0  30 ‘ ’ 

3  4  3  8 

’  3’6 

CV5HOCt=AHAX1 ( ( EPPP/£FP2B-PPpSAV/CD2PSAVI , 

0  4  ?  9 

*  (  £F.CP/P=C  0  P-P"PSAV/PC2PSA  V  )  ) 

0440 

:  =  (CVS-C.1T  GVEAO.'GO  *0  3  ’  05 

044  1 

‘  3  1 1 

NS'.'B  =  «AX0  ;  2  .  (  IP!  v;  50  .  »OVSHOOT  )-49)  ) 

0442 

GO  T0  3107 

0443 

’  3 ' 9 

CONTINUE 

044  4 

C 

0445 

c 

COLLAPSE  YIELDING  ONLY  STRESS  INCREMENT 

0446 

c 

0447 

'600 

AWC*AWC*DWC 

0446 

CALL  MAD0(ARAEC.6,1, DGCDS , 1 .O.DLAMC.ARAEC) 

0449 

CALL  MMULT(ACE. S . 6.  DGCDS. 6 .1.1.0. ACE3 ) 

0450 

CALL  MMUlT( ACEG. 6 .1.81,1,6,1.0, ACEP) 

0451 

CALL  MAOO( ACE . 6 , 6 , ACEP , 1  . 0 . -  1 . 0  ,  ACEP ) 

0452 

CALL  MMULT (ACEP . 5 , 6 , DEP , 6 . 1, 1 .  O.DSIG) 

0453 

CALL  MADO(SIG . 6 , 1 , DS I G . 1 . . 1 .  ,  S I G > 

0454 

WTYPE-1 

0455 

GO  TO  2000 

0456 

c 

103 


»  .'-  .vV' 

►  ■  M  /k. 


EXCLUSIVE  VIELOING  ONLY  STR£S3  INCREMENT 

AWP  =  AWP  +  DWP 

CALL  V40C(ara£P,6, 1 , OGpDS , 1 .0.CL4MP.ARAEP) 

CALL  Y»U  L T  (.  ACE  ,6,6,  DGPDS ,  6 , 1  ,  1  .  3  .  ACEG  ) 

CA.L  wv J LT ( AC EG,6,1,B2.1.6,’.C.ALEC) 
call  MA00(ACE,6, 6.ACEP, 1 .0,-1 .O.ACEP) 

CALL  «*"ULTC  AC  EC,  6 , 6  ,  DEP ,  6  ,  1  .  1  . C , OS  I G  ) 

CALL  MA DD(  S I G .  6  .  1  .  DS I G  ,  '  .  ,  '..3:5.' 

MTVpg-J 

-c  2Dcc 

C 0M9 I  NED  COLLAPSE  AND  EXPANSIVE  YIELDING  STRESS  INCREMENT 


4  r  i 

i  ?  2  2  AwC  =  A«*C ♦ 

DaC 

A  7  2 

r-wc- 

•  ■  } 

4  V4" 

-  ;  4 i a  sr  , 

•  -•  ;  '  *  : 

*.  V  *  L  ~  & 

i  ■*  4 

.  -  _  L  v  4  7: 

^ ;  i,E'4[P( 

c. 

#  z  3  :  7: 3  ,  * 

'  .  “•  _  i  **  ►  ,  A  4  4 

*  -  5 

j  _  _  V  V*  ■ 

. ' . AVE,  5 

c 

.  :?a^ ,9.2. 

■  : . -cep; 

4  **5 

'  A  w  L  v  v  j 

_'r  '  -C  , 

5  . 

2,95t^,2,f 

: ,  a  c  e  ^  2 

;  -  - 

v  4 

c ;  - :  s .  e , 

a  c  :  c ,  - 

4  “*  -9 

.  A  L  •_  vVl 

_  "  :  A  C  E  P  . 

5  ,  C  £  R  .  5  ,  ■  , 

•  “  .75’  7o 

A79 

Call  va; 

D; 3 !G. 6 , 

03 IG,  i  ,  ’  .  , 

,  SCO; 

~  E" E  Pv I N£  EnC-CP  -  I  N.  R£V£N ’■  .  3  .  S  : \C  PE«En  '  ;■  375£3S  S*i*E 

all  a  l  t ;  *i  / <_  ■- ;  g  .  s :  e  v .  s "  * c ,  ; ,  : 

•>ECK  FOR  PR  CONSISTENCY  -  Sl'BC-'CLE  IP  NOT  ADEQ.'ATE.Y  CA’ISEIEO 

:  A '.MTvPE  . .  2  AND.  NSL'8.E-5-C)G3  "C  SCC0 

; p(NS.B . EC . 0 )GO  TO  3100 

:Sw-e»isuB*i 

I  c  (  I  S'jS  .  GT  .  NSJ6  )GO  T0  5000 
CUECK  SPA  -L  FOR  EACH  SL‘9 -  INCREMENT 
I  z  i  SOCT .GE.-(APEX) )GO  "0  3015 
GO  T0  6001 
GO  -0  TOO 

/  I E  *  E  I N  E  NO.  OP  SjBCYCLES  BASED  CN  RATE  OF  C-ANGE  Oc  O'-iPun 
gepPccpp^sOCT, TOCT . COS  3 W ) 

£rc;.p=pp2F  ;  4«p) 

:=;AImO.GT  A^OPY )BFP2P=BFP2P» AWSOFt»FP2P( AWPP* ) * ( ; . C-AWSGFT ) 

ppptcc-ccp2P-9FPP 

IFfFPDIFF.GE. -(TCL2) )GC  TO  5000 

AwOSAV-ST.ATSA  V<  3  ) 

0««'*0PP2PDW(AWPSAv) 

C*.7  »0pb2PDWl  awc  ) 

I  F(  0*2  .  GT  .  S«lNU»*;GO  TO  3104 
I  F  (  DW  1  .GT  .SMLNUM)N SUB-1  0  00 
I F ( DW 1 . LT . SMLNUM jNSLB«2 
GO  TO  3109 
OV3HOOT«OW 1/DW2 
I F( OVSHOOT . LT . 640  .  )GO  TO  3106 

Nsue» i ooo 

GO  to  3  ’  0  9 


>5K5 


A\V>  \ 
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3106  NSU8-IFIX(50 . *OVShOOT)-49 
3  1  C 7  NSU8*MIN0(NSU6  .  ’COC) 

3 ' 08  I F(NSU8 . LE .  1  )G0  TO  5000 
3  0  09  R5U8-1  /NSUB 

CALL  HSCALRv  PEPI  .  6  .  1  .  RSL'8  .  OEP) 
CA.L  “5CALR(SIGI .6. 1 . 1 .C.SIG) 
3300  DO  33:0  1-1 ,NSVAR( 10) 

3  3 ’ 0  STATE; I )*STATSAV( I ) 

isue«i 

AGOC'r  =  3‘rAT  E<  1  ) 

A*C«STA"  ■;  0  ) 

4«P*S1 ATE ( 3 ) 

AJ  2-ST  AT£ ( 4  ) 

ASPALL-S’4*£'5) 

ACC'NPC’  =  5TA'rE  ;  6  ) 

4  5V AP-ST A  *  =  t 5 ) 


■ensile  c4!..;p; 


S'^ESSES  "NO  ’  4  . 


c' 4  =  e ■; )go  To  too 

s cc '  cont:nue 

5  0  0  3  DO  5  3  05  I  *  ’  .  3 

-304  3:G(0j*-454i 

5  3  35  3 ! G ; l ♦ 3 ) * 0 . 0 

:  - :  \ s .  5  i»t  .  0  ;gc  *  o  f  1  ? 5 

as  pa  „_«A3=e4v»  •  ;  *de  5 :  o  1  -  zi~ :  *  3 ) 

*  c  c  *  .*  7 

-33:  A  3  A  A  L  t_  -  4  3  -  E  a  •/  »  4  » ;  %  i  3  .  -3  .  ;  N  3  .  r  -  I  3- .  8  ♦  ’  ,  •  O'  E  ' 


EnERGv  CEPE '<  DEN*  PRESSURE  C  0Nt  s  3  9  j  t  I  C'N 

’uE  34v  pel. at : 3ns-“ i p  :s  e*p:r;c4.  w:->  vn ; T - de^en o^en *  :ons*an*s 

UNITS  ARE  GMS.  TERGS.  vgAR  (C-S-mS)  .  DENSITY  anC  E  E  5  3-  *  A5E 
CON  V  E  "  E  D  TO  fv  /CO  AND  Tesgj.-JM  £>o;q9  *  n  EN  *  Ec  I  NG  tp  3  :  ECU  A  T  3  ON 
GAM  :*3ELC  IS  01 VENSICMLESS . 

’333  IF( A«£ . N£ . ’ . 3 )G C  TO  ’900 
?  3  5  0  DpSOL  E  .  c  0 
GPv 45R- j  ,  C 
OCvapEtO  0 

:c(EcAC.GE.O.C )G0  TO  7300 
eex«exp(epa:> 

E5TAP»'vENGN-AE“ElT  ;*(  1  .  -  £  EX  ) 

00 « AM ax  1 (0“N . OENS*OENCONV(  I  UN  I TS) ) 

EE*AM4X 1 ( ES’AR. A£MELT)/( 00*ENGC0N V ( IUNI TS) ) 

Gam. ( o  3  5  *  A  _0G 1  C ( E  6 ) -0  464)*»2  ♦  0  AO  *  0 .  1 2 • A LOG ’ 0  (  DO ) 

»5 vAP»GAM»06NS*ESt AR 
7  3C0  CONTINUE 

MODIFY  STRESS  TENSOR  TO  REFLECT  mo  EOS: 

(’)  ADO  VAPOR  PRESSURE  TO  STRESS  TENSOR 

(2)  REDUCE  DEVIATOR  STRESSES  OUE  TO  PARTIAL  OR  FULL  MELTING 
AND  RECALCULATE  STRESS  TENSOR 


7A00  00  7420  1-1.6 


g 

1 

057  ’ 

0572 

0573 

0574 

ft 

7 

OS’S 
0576 
0571 
05  78 

V 

0579 

rj> 

«  C  O  rt 
kwOW 

058-> 

s:ev( : ;«soev( i )*4»axi (e*ac . o . c : 
s:g( i ;«sdev( i )»(soct+apvap)*amcol(  i  ) 

8*ODN»(SOCT+APVAP-SOCTSAV)/OEV 


7420  5IG(i;«S 
8*ODN»(S 

7900  CONTINUE 


CONTINUE 

S.  **  ST  BA  I N  COM°ON£NTS 

C  A L  HAD0(APA  =  T,5,  1.0EP.  ’.0.  1  .  C  .  A  R  A  £  T  ) 
CALL  *AOD-;  APAEt.  6  .  1  .  ARAEC  .  1  !  0  .  -  ’•  .  C  .  ABA  EE  ) 
Call  I*ADD(  AAAEE  ,  6  ,  1  .  ARAEP .  I  .  0  .  '  1  C.abaeE) 

cont:nue 

Save  maximum  cl  spent  vosllii  c's  col san’ 
TI^ES'EB  ;a_CJ.4*!Cn 

:  -  .  v«  -  y  _•  t  _  z  ■  _';30  "  QCC- 

*  '*  I  r  *  A  m  a  ■  a  c  i  ^ j  A  ^  £  z>  t  j  •  ' 

...  v  1  N  -  A  m  l  >  '  ;  7  v  ~  ^  \  .  5  C  E  -  '  3  3  '  ; 

' c  ssec 

CVCCA*4CE(  ■  .  1  ; 

eM;:DN  =  AM4x  •  (  A^C  .  (  C^CCN*  v  1  -  A  PCI  )/;  3*  (  :  -  AD 

•'  vr.?A'sAM4k  *  *+ 2  ,  K  CVCDN  *  '  1  J  ^  pfj  *  •  f  t  2*  (  '  - 

“‘‘A 

-1--  V2C  A.  S  :  5 1 G  L  4  )  .  :■ .  •  .  .  "  g  P  '  ,  i  I  ' .  4  )  ,• 

CA..  vsca..s;s:c,s.  •:  . s i c ; 

-  A  -  w  l':Ci.s.:-3!.5.',.g:,;E5;j 

SAVE  S'A’E  VARIABLES 
S *  a T E  (  1 >  =  ASOCT 

sT a ' e ( :  )»a*c 
S*ATE(  3  ;  =  A«P 
ST  AT £ ( 4 ) ■ AJ  2 
S’A’EC5)*A3pall 
S*A"ECS;  -Jw>;ee 
S*a'e;7)*aev 
S'A^e  ;  5  I=ASr’ENG 
S’ATE( 3)»APVAP 

RETURN 


Xversion  semi ink 

updates  used  by  applied  research  assoc.,  albuquerque,  nm 
*1  none* 

Kn  none* 

K/ 

X y  this  file  contains  the  atlas  update  cards 

* /  for  stealth2d  version  4-ld,  as  specified: 

*/ 

«/  bigxxx  -  grid  size  adjustments  to  common  blocks 

*y  addxxx  -  add  variables  to  stealth  common  blocks 

X/  semxxx  -  stealth  2d  link  with  soil  element  model 

*y  semlnk  -  programmatic  interface 

xy  seminp  -  mods  to  material  input  cards  for  sem  input 

xy  semstn  -  engineering  strain  increments 

*/  semmod  -  zonusr  routine  for  calling  sem  models 

x/  semovb  -  overburden  initialization  using  sem  models 

X/  gamxxx  -  substitution  of  gamma-law  gas  for  jwl  eos 

X/  bswxxx  -  boundary  control  switches  -  problem  dependent 

x/  prpxxx  -  pressure  loading  function  in  'myfno' 

x/  lysxxx  -  lysmer  non-reflecting  boundary 

X/  actxxx  -  activity  check  based  on  volumetric  strain  rate 

*y  matxxx  -  matinp  changes  for  resart  material  model  changes 

X/  sbexxx  -  error  flags  for  subcycling  control 

X/  sldxxx  -  corrections  to  slideline  search  algorithms 

*/ 

xy  file  created  :  13  december  1984 

xy  last  changes  t  30  july  1985 

xy 

xy  use  oldps  and  binary  for  stealth  version  4-ld 

xy 

xy  stealth  version  4-ld  is  located  on  the  following  cfs  path 

xy  y0000536ytwodystealth/ccyseismic2 .xx 

xy 

xy  contact  bill  dass  for  info  on  this  file  (919)876-0018 


xy  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xy  xx  xxx 

xy  xx  file:  bigxxx  xxx 

xy  xx  xxx 

xy  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xy 
xy 

xy  adjustment  of  maximum  grid  size 

xy 

»d  maxiyl 
65 

xy  51 
*d  maxjyl 
45 

xy  4 

*/ 

xy 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXKXXXXXX 
xy  xx  xxx 

xy  xx  file:  addxxx  xxx 

X/  XX  xxx 

xy  XXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/ 


X/  update  stealth  zonary/var  to  include  more  extra  variables 

X/  and  to  track  original  zone  dimensions 

X/ 

Xi  zonary2/16 

K  ex7  C&maxiS ,  4) ,  ex8 ( 8maxi& ,  4 ) ,  ex9(&maxi&,  4)  , 

X  xd0(&maxi&,4) ,  ydO ( &maxi&,  4 ) , 

*i  zonvar/13 

X  ex7o> ex7n, ex8o, ex8n, ex9o, ex9n, xd0o,xd0n,yd0o,yd0n, 

«i  gptnew2/117 
c  ex7(ic, jb)=0 . 0 

c  ex8(ic,jb)=0.0 

ex9(ic, jb)=0 . 0 
Xi  xonold2/78 

ex7o=ex7(ic#ic) 
ex8o=ex8(ic, jc) 
ex9o=ex9(ic, jc) 
xdOo=xdO(ic, jc) 
ydOo=ydO(ic,  jc) 

Xi  zonnew2/98 

ex7n=ex7o 
ex8n=ex8o 
ex9n=ex9o 
ex7(ic, jb)=ex7n 
ex8(ic, jb)=ex8n 
ex9(ic, jb)=ex9n 
xdOCic, jb)=xdOo 
ydOfic, jb)=ydOo 
X/ 

X/  update  matary  to  re-include  heat  variables 
x/ 

Xi  matary/21 

common/matary/ 

X  macon(lO), 

X  acon0(10),aconl(10),acon2(10),acon3(10),acon4(10), 

X  acon5(10),acon6(10) , acon7 (10),acon8(10) , acon9( 10 ) , 

x  mashc(lO), 

X  ashc0( 10 ) , ashcl ( 10),ashc2( 1 0 ) ,ashc3( 10 ) ,  ashc4(  10), 

x  a she 5 ( 10), ashc6 ( 10 ) ,ashc7 ( 10 ) , ashc8 (10) ,ashc9( 10), 

x  maexr(lO) 

Xi  matvar/27/ 

X  mexr, 

X/ 

X/  KXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/  XX  XXX 

x/  xx  file:  semlnk  *xx 

X/  XX  XXX 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXX 


x/  xxxxx  afwl  stealth  version  4-ld  link  with  soil  element  model 

x/ 

X/  soil  element  model  common  blocks 

x/ 

Xstring  constant 

common  /constant/wonthd, toothd,sqrt2, sqrt3,sqrt6 ,sqrt23» 


Xstring  number 

common  /number/ 

1 

2 


twopi , bigpos, bigneg 

itest , iundr, iskip, j print , jplot, iplt( 30 ) ,  kount, 
iexec,iprob,nprob,isave(10), jsave, nparamt 1 1 ) , 
nsvar( 11 ) , patm(8 ) , nf ini , nfin2, nf in3, nf in4, nfin5. 


rw 
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3 

Kstring  energy 

common  /energy/ 
Kstring  tens 

common  /stress/ 
common  /strain/ 
Kstring  strinc 

common  /strinc/ 
Kstring  allmod 

common  /allmod/ 

1 

Kstring  propO 

common  /propO/ 
Kstring  propl 

common  /propl/ 
Kstring  prop2 

common  /prop2/ 
Kstring  prop3 

common  /prop3/ 
Kstring  prop4 

common  /prop4/ 
Kstring  prop5 

common  /prop5/ 

1 

Kstring  prop6 

common  /prop6/ 

1 

2 

3 

Kstring  prop7 

common  /prop7/ 

1 

Kstring  prop8 

common  /prop8/ 

1 

Kstring  prop9 

common  /prop9/ 

1 

Kstring  proplO 

common  /proplO/ 

1 

2 

3 

4 

Kstring  proplla 

common  /propl 1/ 

K/ 

X/  update  stealth 
xd  stealth/1-6 


nfotl,nfot2,nfot3,nfot4,nfot5 


engd 

sigx, sigy , sigz, sigxy, sigyz,sigxz, p, ep, s( 3) 
epsx, epsy, epsz , epsxy , epsyz, epsxz , eC  3) 

depsx, depsy, depsz, depsxy , depsyz, depsxz 

ieos, iuni ts>  atmo , conv , rhoref » mtype, gr , grp, 
bsp( 3) , vstnC  3 ) , bmodn, gmodn, state ( 10) 

mloc( 11 ) , prop( 1 ) 

bulk3,shear3 

bulkl , bulk2, cl , shear 1 , shear2, csl 

bulk, shear , ca ,cb,cc,cam,tcutl,fcutl, rule, esp 

rub, ru, game, tauc, tauy, gmax, gamo, bulk4 

aki,akl,ak2,akim,aklm,ak2m,agi,agl,ag2,ag3,ag4, 
ac , am, bb, ccc , a  ri , a r 1 , ar2, aw, ad, akh,akhm, el 

rnls, rnus, bkl (8) , ebl (8 ) , pol (8 ) , bku(8 ) , 
pbu( 8 ),pou(8),stl,yl,sl,vml, fstype, st2, y2, 
s2, vm2,sptype,pxcut, pycut, pzeut, rihe, emelt, cte, 
pbl (8 ) , ebu(8) , evmax, evgrav , espal 1 , epl ( 3) 

ekur, en, pois, c, pc, etal , curvm, r ,ss, t, alpha , 
beta , pw, elw,q, wppk,a,b,sigma3,wc,wp 

hk,hkur,bn,hc,hrf,hphi,hkb,hkbur,hg,hf  ,hd, 
hsigma3, hei , heur , hfsdif f , hnui , heamax, hevmax 

exa , exw, ejrl , exb, ej r2 , vdet , rlvc j , 
tact, dact, exe, exv 

akur,an,apoi ,acrv,aacc( 4) ,aapc( 4) ,abrx(3) , 
ar ,aey,amy,aetal ,apbar ,al ,aalph,abeta, 
a t g, a r g, as g, apex, ahswtch, a ha, a hn, a hi am, a hgam, 
ahbet ,aconfp,asoct,acc,apc, awppk  >aq,aa,ab, 
aeta20,awc,awp,fp2ptyp,awsoft,aj2,aetar,ahy 

pcon( 175) , wgmax, emuz, umax, umin, gpsz, vole, eng 

program  card  to  include  sem  i/o  units 


program  stealth(stlinp, tape5=stl inp, stl out , tape6 =stl out , 

K  fiche,tape20=fiche,stdlib,tape31=stdlib, 

K  usrlib,tape32=usrlib,matlib,tape37=matlib, 

K  tape7 , tape8 , tape9 , tapel 0 , tapel 1 , tapel2, 

K  tapel3,tapel4,tapel5,tapel6,tapel7) 


c 

Ki  stealth/73 


insert  sem  common  into  stealth  as  contiguous  block 


tnumberS 
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Aconstant & 

Aprop0& 

AproplS 

Aprop2& 

Aprop38 

Aprop4S 

Aprop5& 

Aprop6& 

Aprop7& 

4prop8& 

Aprop9S 

AproplOA 

ApropllaS 

c 

c 

c  initialize  soil  element  model  constants  and  unit  numbers 
c 

wonthd=l ./3 . 
toothd=2 ./3 . 
sqrt2=sqrt(2 . ) 
sqrt3=sqrt( 3 . ) 
sqrt6=sqrt(6 . ) 
sqrt23=sqrt( toothd) 
two pi =8 . XatanC 1 . ) 
bigpos=l . 0e90 
bigneg=~l . 0e90 
c 

nfinl  =  5 
nfotl  =  6 
nfin2  =  37 
nfin3  =  38 
nfin4  =  39 
nfin5  =  40 
nfo±2  =  6 
nf ot3  =  39 
nfot4  =  41 
nfot5  =  42 
iexec  =  1 
itest  =  10 
c 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/  XX  XXX 

X/  xx  file:  seminp  xxx 

X/  XX  xxx 

X/  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X/ 

X/ 

X/  input  changes  for  sem  material  models 

x/ 

xi  matset/15 
c 

Azonary2A 

AnumberA 

AallmodS 

SpropOS 

c 

*d  matset/19 

X  namcap(80 ) ,  namengC80 ) , namsem( 80 , 11 ) 

*d  matset/51-62 

c  set  up  names  for  sem  models  mmdl=-l  thru  -11 
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data  (namsem(i , 1 ) , i  =  l ,80 )/ 

*  lhe,  lhl,  lha,  lhs,  lht,  lhi,  Ihc,  73*lh  / 
data  (namsemCi, 2), i=l, 80)/ 

*  lhv,  lhi,  lhs,  lhc,  lho,  75*lh  / 
data  (namsem(i,3),i=l,80)/ 

X  lhe,  lhl,  lhp,  lhl,  lha,  75Xlh  / 
data  (namsemCi, 4), i=l, 80)/ 

X  lhp,  lhy,  lhk,  lhe,  76*lh  / 
data  (namsemCi, 5), i=l, 80)/ 

*  lhc,  lha,  lhp,  77*lh  / 
data  (namsemCi ,61 , i=l , 80 )/ 

*  lha,  lhf ,  lhw,  lhl,  76*lh  / 
data  (namsem(i,7),i=l,80)/ 

X  lhl,  lha,  lhd,  lhe,  76*lh  / 
data  (namsem(i , 8 ) , i=l, 80 )/ 

x  lhh,  lhy,  lhp,  lhe,  lhr,  75*lh  / 
data  (namsemCi, 9), i=l, 80)/ 

*  lhj ,  lhw,  lhl,  77xih  / 
data  (namsemCi , 10) , i=l , 80 )/ 

*  lha,  lhr,  lha,  lhl,  76»lh  / 
data  (namsemCi , 11 ), i=l , 80 )/ 

x  lhw,  lha,  lhg,  lhn,  lhe,  lhr,  74*lh  / 

*d  matset/143-210 
c 

2000  imdl  =iabs(mmdl) 
ieos=imdl 
maeos(l)=9 
mayld( 1 ) =9 
mashr(l)=9 
mmmatC 1 )=-l 
mmat=mmmat( 1 ) 

2100  do  2105  ln=l,80 
2105  nammdl(ln,l)  =  namsemC In, ieos) 
c 

Xd  matset/224 
3100  call  inpeos 

read(nfinl , 3105) (nammat( 11 , lmat), 11=1,80) 

3105  format(80al) 
c 

c  store  reference  density  and  initial  geostatic  pressure 

c 

ardn( lmat)=rhoref 
a c oh ( lmat, 1 )=grp 
c 

c  store  model  parameters  (80max)  in  matary 
c 

3200  do  3205  ii=l,nparam(ieos) 
rskip=(ii-l)/10 
iiskip=ifix( rskip) 
ifCiiskip.ge . 4)iiskip=iiski p+1 
locary=10x(ii-l+iiskip)+lmat 
3205  aeos0( locary)=prop(mloc( ieos)+ii ) 
c 

X/ 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
K/  XX  XXX 

X/  **  filet  semstn  *** 

x/  xx  xxx 

X/  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
Xi  zonstn2/66 


W.v 


SstrincS 

Xi  zonstn2/108 

c 

c  calculate  engineering  x-y  strain  increments  for  sen  models 
c 

facx=0 . 5*dlth/xd0o 
facy=0 . 5*dlth/yd0o 

depsx=facxX( (xvlhbi — xvlhtl )  +  (xvlhti — xvlhbl ) ) 
depsy=facyX( ( yvlhtl-yvlhbr )+(yvlhtr-yvlhbl ) ) 
depsxy=facx*( (yvlhbi — yvlhtl )  +  (yvlhti — yvlhbl ) )  + 

X  facy*( ( xvlhtl -xv lhbr )+Cxvlhtr-xvlhbl ) ) 

depsyz=0 . 0 
depsxz=0 . 0 
xi  zonstn2/139 

c  translational  plane  strain 
depsz=0 . 0 
xi  zonstn2/143 
c  axial 

depsz=vsrhXdlth-depsx-depsy 

c 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/  XX  XXX 

x/  XX  file:  semmod  xxx 

X/  XX  xxx 

X/  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xi  zonusr/15 
c 

c  make  call  to  sem  model  in  zonusr 

c 

4  matary  8 

StimvarS 

SnumberS 

4tens8 

Sstrinc8 

4energy& 

Sal lmodS 

SpropOS 

c 

c  array  locations  for  extra  (state)  variables  exl-9  in  zonvar 
dimension  exg(18) 
equivalence  (exg( 1 ) , exlo) 
c 

xd  zonusr/22-41 
c 

ieos=lmdl 
rhoref =rdn 
c 

c  set  material  properties  in  /prop/  from  /matary/ 
c 

200  do  210  ii=l >nparam(ieos) 
rskip=(ii-l )/10 
iiskip=ifix( rskip) 
if ( iiskip . ge .  4)iiskip  =  iiskip+l 
locary=10X(ii-l+iiskip)+mpno 
210  pr op ( ml oc(ieos)+ii ) =aeosOC locary) 
c 

c  initialize  material  model  state  variables 
c 

220  do  225  ii=l>nsvar(ieos) 

jj=iix2-l 


-v.% -»\v.v  *.  .%■ 


-V  N  .VA  . 


c  check  for  maxsbc  exceeding  allowable  limit 

260  if (maxsbc . It . 30 )  go  to  270 
write(nfmsg, 1261)  maxsbc 
writeCnfprt, 1261 )  maxsbc 
1261  format 

x  (10x,10hxxxxxxxxxx 

x  /10x»10h  error 

X  /10x,27hmax  subcycle  limit  exceeded 

X  /10x,8hmaxsbc=  ,ilC) 

nserr  =  4 
nsext  =  4 
go  to  990 

270  dltn  =  dltsbc 

X/ 

X/ 

X/  KXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/  XX  XXX 

x/  XX  file)  sldxxx  xxx 

X/  XX  xxx 

X/  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X/ 

X/ 

X/ 

x/  slider  modifications  83/03/16/14:05 

X/ 

X/  debug  of  minimum  search  routine 

x/ 

xd  fndpta2/42 

if ( lpta . ge. 1 .and. lpta . le. nlopp)  go  to  50 
c 

c  do  minimum  distance  search  of  entire  opposite  grid 

dstsml  =  0.5  X  bignum 
do  25  11=1. nlopp 

dsttry  =  (xpnf  -  xpnargf 1 1 . i opp) )  XX  2 

X  +  (ypnf  -  ypnargC 11 , iopp) )  XX  2 

if (dsttry. gt. dstsml)  go  to  25 
dstsml  =  dsttry 

lanow  =  11 

25  continue 
go  to  400 
50  continue 
Xi  supget2/45 
c 

c  set  left  and  right  points  to  current  point 


c 

go  to  30 
c 

c  normal  loop 
c 

27  continue 


c 

x/ 


K/  XXXXXXXXXXXXKXXXXXXXXXXXXXXX.  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKX 
X/  XX  XXX 

*/  **  filei  matxxx  xxx 

X/  XX  xxx 

X/  XXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

XX 

Xd  matinp/98-99 

mmmdl ( 1 ) =  mmdl 
if (mmdl . eq . 0)  go  to  90 


atlas  input  for 


modifications  to  matinp  to 

allow  changing  of  material  models  during  restart. 


file  last  changed  83/06/05/13:58:36 


x/ 

X/ 


update  for  minimum  subcycle  time  step  and 

maximum  subcycles  checks  and  stops 


file  last  changed  83/03/15/20:20 


X/  XXXXXXXXXXXXKXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXKXX 
X/  XX  xxx 

*/  **  file:  sbcxxx  xxx 

X/  XX  xxx 

X/  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

X/ 

Xd  sbctim2/79 

200  if(dltn.ge.dltsbc)  go  to  270 
»d  sbctim2/83 
c 

check  for  subcycle  time  step  less  than  minimum  time  step 
if(dltsbc.ge. dltmin)  go  to  260 
write(nfmsg, 1251 )  dltmin,  dltsb- 
writeCnfprt, 1251 )  dltmin,  dltsbc 
1251  format 

( lOx, lOhxxxxxxxxxx 
/10x,10h  error 
/10x,23hsubcycle  time  step  less 
/10x,19hthan  allowable 
/10x,8hdltmin=  ,lpel2.5 
/10x. 8hdltsbc=  ,lpel2.5) 

=  9 
-  9 


250 


x 

X 

X 

X 

X 

X 


nserr 

nsext 


go  to  990 
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fl 

£ 

I 


V, 

i 


R 


II 

> 


c 

c 

c 


c 


rlv(ic,jb)  - 
vsr(ic,jb)  = 
com(icob)  = 
dns(ic,jb)  = 
zms(ic,jb)  = 
zie(ic,jb)  = 
zde(ic,jb)  = 
zse(ic,jb)  = 
zke(ic.jb)  = 
prh(ic,jb)  - 
avs(ic,jb)  = 
sss(ic,jb)  = 
trv(ic,jb)  = 
txx(ic.jb)  = 
tyyCic,  jb)  = 
tzz(ic,jb)  = 
txy(ic,jb)  = 
sxx(ic,jb)  = 
syy(ic,jb)  = 
szz(ic,jb)  = 
pxx(ic,jb)  ~ 
pyy(ic,jb)  = 
pzz(ic,jb)  = 
pxy(ic,jb)  = 
epi(ic,jb)  = 
eps(ic,jb)  = 
yld(ic.jb)  = 
shr(ic,jb)  = 
indCic, jb)  = 
mpn(ic,jb)  = 
ign(ic,jb)  = 
bf s ( ic,jb)  = 
act(ic,jb)  = 
qda ( i c ,  j  b )  = 
txl(icjb)  s 
ex2(ic,jb)  = 
ex3(ic,jb)  = 
ex4(ic,jb)  = 
ex5(ic,jb)  = 
ex6(ic,jb)  = 
ex7(ic,jb)  = 
ex8(ic,jb)  = 
ex9(ic,jb)  = 
xdO(ic,jb)  = 
ydO(ic,jb)  = 

energy  check 


rlvCic, jc) 
vsrC ic, jc ) 
comCic, jc) 
dnsCic, jc) 
zns(icoc) 
zieCic, jc) 
zdeCic, jc) 
zseCic, jc) 
zkeCic, jc) 
prhCic, jc) 
avsCic, jc) 
sssCic, jc) 
trvCie, jc) 
txx( ic, j  c ) 
tyyCic, jc) 
tzzCic, jc) 
txyCic, jc) 
sxx(ic> jc) 
syyCic, jc) 
szzCic, jc) 
PxxCic, jc) 
pyy(ic, jc) 
pzzCic, jc) 
pxyCic, jc) 
epiCic, jc) 
epsCic, jc ) 
yld( ic  >  jc ) 
shr( ic, jc ) 
ind( ic,  jc) 
mpn( ic, jc ) 
ign(ic, jc) 
bfs(ic, jc) 
■  ct( ic, jc  ) 
qdaCic, jc) 
exlCic, jc) 
ex2 ( i c , j  c ) 
ex3(ic, jc) 
ex4(ic, jc) 
ex5(ic, jc) 
ex6 ( ic, jc ) 
ex7 ( i c , j  c  ) 
ex8(ic, jc) 
ex9(ic, jc) 
xdOCic, jc) 
ydO ( ic , j  c ) 


fac  =  zms(ic,jb)  /  rdn 
tie  =  tie  +  fac  X  zieCic,jb) 
tde  =  tde  +fac  *  zde(ic,jb) 
tke  =  tke  +  zkeCic, jb) 
tms  =  tms  +  zms(ic,jb) 

txm  -  txm  +  zms(ic,jb)  *  (xvl(ic.jb)  +  xvl(il,jb)  + 

*  xvl(ic,jw)  +  xvl C i 1 , jw) ) 
tym  =  tym  +  zms(ic,jb)  *  Cyvl(ic,jb)  +  yvl(il,jb)  + 

*  yvl(ic,jw)  +  yvl(il,jw)) 
teg  =  tie  +  tke 


call  timpro 
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WV3V 


.\s\ 


e 

*i  propro2/58 
c 

tnatvar& 

c 

igptary2& 

c 

*i  propro2/182 
c 

c  sava  the  old  activity  check 

c 

noro=nor(ic, jc) 
if (noro . gt .  4 )actblc=actsav 
if (noro .gt .  3)actsav=act(ic> jb) 
c 

c  only  test  interior  points 

c 

if Cnccyc. le. 3)go  to  27 
if (noro.ne.5)go  to  27 


test  for  activity 


actsum  = 
if  + 

K  + 


act( il , jt) 
act(il« jc) 
actblc 


if (abs (actsum) . gt . 1 


+  act(ic,jt)  +  act(ir,jt) 
+  act(ic,jc)  ♦  act(ir,jc) 
+  act(ic,jb)  +  act(ir,jb) 
,  Oe-15)  go  to  27 


no  activity  -  set  variables 


/gptvar/ 


c8027 


write(6 ,8027)nccyc, i , j ,  actsum 

f ormat( 'act . . . nccyc= ' » i3»  '  i,j',2i3,’  actsum5 * , IpelO . 3) 


xpnotl 

ypnotl 

xpnobl 

ypnobl 

xpnotr 

ypnotr 

xpnobr 

ypnobr 


xpnotr 
ypnotr 
xpnobr 
ypnobr 
xpntic, jc) 
ypnfic, jc) 
xpndCf  jb) 
ypn(ic, jb) 


set  variables  in  /gptary/ 


xpn(ic, jb) 
ypnfic, jb) 
zpn(ic, jb) 
xvl(ic,jb) 
yvl(ic, jb) 
zvlCic,jb) 
xac(ic.jb) 
yacCic, jb) 
zac(ic, jb) 
gxk  C  i  c ,  j  b ) 
gyk ( i c , j  b ) 
nor( ic, jb) 
nbm(ic» jb) 
nbv(ic> jb) 


xpn(ic>  jc) 
ypn( ic»  jc) 
zpn( ic>  jc) 
xvl(ic»jc) 
yvl ( ic, jc) 
zvl(ic,jc) 
0.0 
0.0 
0.0 

gxk  ( i  c ,  j  c ) 
gyk  C  i  c ,  j  c ) 
nor( i c , j  c ) 
nbm( ic ,  j c ) 
nbv(ic, jc) 


set  variables  in  /zonary/ 
dll(ic>jb)  =  dll(ic,jc) 


A  >1, 


|L  it 


delx-0 . 5#(xpn( iz, jz)-xpn( iz-1 , jz)+xpn(iz, jzb)-xpn(iz-l ,  3 zb) ) 
dely=0 . 5x(ypnCiz, jz)-ypn( iz-1 , jz)+ypn(iz»  jzb)-ypn( iz-1 , jzb) ) 

gat  sound  spaads 

sssprh  =  sqrt(sss( iz/ jz) ) 
ifU.eq.l)  rdn=1750. 
rhoref =1900 . 

sssshr  =  sqrt(shr(iz, jz)/rhoref ) 

set  x  and  y  direction  dll  and  ssp 

dllx=abs(delx) 

dlly=abs(delx) 

sspx=sssprh 

sspy=sssshr 

compute  lysmer  boundary  terms 

200  trmxl  =  sspx  *  dnsb  /  dnso  *  dlto  /  dllx 
trmyl  =  sspy  *  dnsb  /  dnso  *  dlto  /  dlly 
trmx2  =  1.0-trmxl 
trmy2  =  1.0-trmyl 
trmx3  =  1.0+trmxl 
trmy3  =  1.0+trmyl 


compute  motion 

300  xvlh=(xvlm#trmx2+xaco*dlto)/trmx3 
yvlh=(yvlmXtrmy2+yaeoXdl to)/trmy3 
if (absCxvlh) . lt.xvlmin)  xvlh=0.0 
if (abs(yvlh) . lt.yvlmin)  yvlh=0.0 
if ( j . eq . ncrow)  yvlh=0.0 

recompute  accelerations 

xaco=(xvlh-xvlm)/  dlto 
yaco  =  (yvlh  -  yvlm  )  /  dlto 

lysmer  boundary  diagnostics 

write(6  »  500 )i , 3 , sspx, sspy.xaco, yaco, xvlh,yvlh 
500  format( 'mybdy . . . i , j ’ , 2i3, *  sspx, y ' , Ip2el0 . 3, ’  x, yaco ' , lp2el 0 . 3, 
500+  •  x,yvlh' , Ip2el0 . 3) 

f  XXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
'  XX  XXX 

'  xx  filet  actxxx  XXX 

'  XX  XXX 

'  XKXXXXKXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


activity  check  based  on  volumetric  strain  rate 
zonstn2/130 

set  activity  switch  in  zonstn 
actn=0 . 0 

ifCabs(vsrh) .gt.l . 0e-15)actn=vsrh 


Xd  myf no/26 -27 

go  to  (100,200)  iftyp 
c 

c  dihest  pressure  function 

c 

100  if(var.gt.fca(3))go  to  130 
c 

120  val=fca( 1 )Xexp(-fca(2)Xvar)x(l . 0-var/fca(3) ) 
go  to  900 

C  130  val=0.0 

go  to  900 
c 

c  speicher-brode  overpressure 
c 

200  call  spbrode(fca(l),fca(2),fca(3),fca(4),nccyc,var,val) 
go  to  900 
c 

900  continue 

X/ 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

»/  V¥ 

x/  xx  file i  lysxxx  xxx 

X/  XX  xxx 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXK 

X/ 

X/ 

X/  lysmer  non-reflecting  boundary 

x/ 

Xd  mybdy2/14 
c 

c  compute  motion  of  lysmer  non-reflecting  boundary 

c 

c 

&matvar& 

c 

IgptvarS 

c 

SzonvarS 

c 

&gptary2S 

c 

&zonary2& 

c 

*timvar& 

Xd  mybdy2/20~21 
c 

c  set  appropriate  zone  indices 

c 

iz=ir 

if(ic.eq.ncgpt)  iz=ic 

32  =  jt 

jzb=jc 

if(noro.ge.7)  jz-jc 
if (noro . ge .7 )  jzb=jb 
c 

c  get  density  and  zone  length 


dnso~dns(iz, jz) 
dnsb=dnso 
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e 

c  dihest  points  may  never  have  negative  x-coordinates 

c 

if(i.gt.l)go  to  108 
if ( j . ge . 5 -and. j . le. 17 )go  to  110 
108  continue 

X/ 

H/  control  of  boundaries  for  pseudo-ld  calculations 

X/ 

X/  *d  gptbdy2/51 
X/  c 

M/  c  a  pressure  boundary  has  been  specified 

X/  c  lnbv  =  1  for  bottom-side  pressure  history  (y-dir  prop) 

X/  c  lnbv  =  4  for  right-side  pressure  history  (x-dir  prop) 

X/  c 

K/  120  continue 

X/  if(lnbv.eq.4)  call  gptprhC val , timo, lnbv) 

X/  if ( lnbv . eq. 4)  ex7n=val 
X/  c 

X/  xi  gptmot2/225 
X/  c 

X/  c  no  x-velocity  at  bottom-right  wall-interacting  grid  point 
X/  c 

X/  if(noro.eq.3)  xvlh=0.0 

X/  c 

X/  Xi  gptmot2/240 
X/  c 

X/  c  no  y-velocity  at  top-left  pressure-boundary  grid  point 
*/  c 

X/  if(noro.eq.7)  yvlh=0.0 

X/  c 

X/  c  no  y-velocity  at  bottom-left  pressure-boundary  grid  point 
x/  c 

x/  if(i.eq.l  .and.  j.eq.l)  yvlh=0.0 

*/  c 
*/  c 

x/  xi  gptmot2/247 
*/  c 

X/  c  pressure  bdy  points  may  never  have  negative  x-coordinates 
*/  c 

X/  c  if(i.eq.l)go  to  110 
X/  c 

*d  gptnew2/102-103 
c 

c  don't  zero  out  zone  interior  variables  used 

c  for  storing  applied  pressure  and  impulse 

c 

ex7 ( ic, jb)=ex7n 

ex8( ic, jb) =ex8(ic, jb)+ex7n*dlth 
c 

*/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/  XX  XXX 

x/  xx  file«  prpxxx  xxx 

X/  XX  XXX 

X/  XXXXXXXXXXXXKXXXXXXXKXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/ 

X/ 

Xi  myfno/21 

data  iftyp  /l/ 

X/ 


I 

K 

8 


& 


S- 

.V 


I 


5? 


*d  zonsss/112-116 
c 

c  sound  speed  squared 
c 

sssprh  =  eosO*(eosO-l)K(zien/dnsn) 
c 

*/ 

x/ 

K/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXX 
K/  XX  XXX 

*/  xx  file!  bswxxx  xxx 

K/  XX  XXX 

*/  XXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

*/ 

X/ 

X/  control  of  boundaries  for  simquake  ii  ff  problem 

x/ 

*i  gptbdy2/25 
c 

Azona ry2& 

SzonvarS 

c 

*d  gptbdy2/51 
c 

a  pressure  boundary  has  been  specified 

lnbv  =  1  for  bottom-side  1  atmosphere  pressure  application 
lnbv  =  2  for  right  side  non-ref lective  boundary 
lnbv  =  5  for  left-side  pressure  history 


c 

c 

c 

c 

c 


120  continue 

if ( lnbv . eq. 1 )  val=101379. 
i f ( lnbv . eq . 2)  val =ex9(i , jx) 
c  i f ( lnbv . eq . 2 )  write(6 > 9000)i , j , val 
c9000  f ormat( ' gptbdy . . . i , j ' , 2i3» •  val= ’ » lpel 0 . 3) 
if ( lnbv . eq. 5)  call  gptprh(val , timo, lnbv) 
if ( lnbv . eq. 5)  ex7n=val 
go  to  121 

121  continue 


s* 

m  d 


Xd  gptmot2/226 
c 

c  no  x-velocity  at  top-left  grid  point 


if ( noro . eq. 7 )  xvlh=0. 
if  (j . eq.45)yvlh=0.0 
if ( 3 . eq. 45)ypnn=ypno 
79  go  to  400 


*i  gptmot2/240 
c 

non-ref lective  boundary  along  right  side 


c 

c 

c 

c 

c 


if ( i . eq . ncgpt)  call  mybdy 
no  x-velocity  at  top-left  grid  point 
if ( noro . eq . 7 )  xvlh=0.0 


■» 

v 


*i  gptmot2/247 


2 
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*  txxn, tyyn, tzzn, prhn, comi ) 

sxxn=txxn+prhn 
syyn=tyyn+prhn 
szzn=tzzn+prhn 
if (grvy . gt . 0 . 0 )ex9n=-txxn 
if (grvx. gt . 0 . 0 )ex9n=-tyyn 

•psn=sqrt(amaxl (atmo/100 . , -2 ./3 . *(sxxnXsyyn+syynXszzn+szzn*sxxn) ) ) ) 
rlvn  =  comi  +  1.0 
9340  do  9350  ii=l,nsvar(ieos) 
j j=ii*2-l 
exg( jj)=state(ii) 

9350  exg( j j+1 )=stateCii) 
rlv(iot)  ~  rlvn 
tyy(i,jt)  =  tyyn 
txxCi,3't)  =  txxn 
tzz(i.jt)  =  tzzn 
prh(i,jt)  =  prhn 
sxxCi.jt)  =  sxxn 
syy(i,jt)  =  syyn 
szz(i,jt)  =  szzn 
exl(i,jt)  =  exln 
ex2(i,jt)  =  ex2n 
ex3(i,jt)  =  ex3n 
ex4(i,jt)  =  ex4n 
ex5(i,jt)  =  ex5n 
ex6(i,jt)  =  ex6n 
ex7 ( i , j  t )  =  ex7  n 
ex8(i,jt)  =  ex8n 
ex9(i,jt)  =  ex9n 
eps(i,jt)  =  epsn 

xd0n=0 . 5*( (xpnntr-xpnnbl )  +  (xpnnbi — xpnntl ) ) 
yd0n  =  0 . 5*( (ypnnti — ypnnbl )+(ypnntl~ypnnbr ) ) 
xdO(i, jt)=xd0n 
yd0(i, jt)=yd0n 

c  if Ci . eq. 2 )write(6 , 8001 )i , j t , ieos,xpnntr, ypnntr ,xd0n,yd0n, 

c  *  trefy, ttempy, grp, txxn, tyyn 

c8001  format( 'i, jt.ieos' ,3i3, *  x,ytr ' , Ip2el0 . 3,  '  x,y0 ’ , Ip2el0 . 3, 
c8001+  / , ' trefy, tty, grp* , Ip3el0 . 3, '  txx.yy' , Ip2el0 . 3) 

Xd  genchk2/357 

ifCnmdl . It. 0)  go  to  380 
Xd  genchk2/378-389 
c 

380  continue 
c 

X/ 

*/ 

X/  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X/  XX  XXX 

X/  **  file:  gamxxx  xxx 

X/  XX  xxx 

*/  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X/ 

X/ 

*/  substitution  of  gamma-law  gas  for  jwl  eos 

X/ 

xd  zonprh/189-191 
c 

c  pressure-dansity-energy  relationship 
c 

prhn  =  ( eosO-1 )Xzien 


X/  XX  XXX 

X/  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X/ 

X/ 

X/  updates  for  overburden  initializationy  using  sem  models 

X/  assuming:  i-lines  are  oriented  in  the  y-direction 

X/  j-lines  are  oriented  in  the  x-direction 

X/  application  of  gravity  stress  is  uniaxial, 

X/  either  in  x  or  y  direction,  but  not  both 

X/ 

X/ 

Xi  genchk2/58 
SnumberS 
Sal lmod& 

ApropOS 
Xi  genchk2/60 

dimension  exg(18) 
equivalence  ( exg( 1 ) , exlo ) 
x/ 

«i  genchk2/103 
c 

c  increment  y-stress  on  row  loop 

c 

trefy  =  trefy  +  ttempyX2.0 
c 

c  initialize  x-stress 

c 

trefx  =  0.0 
ttempx  =  0.0 
c 

Xi  genchk2/l 04 

c 

c  increment  x-stress  on  column  loop 

c 

trefx  =  trefx  +  ttempxX2.0 
c 

Xi  genchk2/262 

c  note:  gravity  stresses  only  calculated  once  per  depth( j-line) 

if ( grvy . gt . 0 . 0  .and.  i.ne.2)go  to  9340 
mnum  -  mpn(i,jt) 

c  assume  soil  properties  for  explosive 

c  if (mnum . eq . 1 )  mnum=2.0 

ieos=-nmndl (mnum) 
rdn  =  ardn(mnum) 
rhoref =rdn 
grp=acoh(mnum,  1 ) 
c  set  matvar  properties 

do  9310  ii=l, 50 
rskip=(ii-l )/10 
iiskip=ifix( rskip) 
if(iiskip.eq.4)iiskip=5 
locary=10x(ii-l+iiskip)+mnum 
9310  prop(mloc(ieos)+ii)=aeosO(locary) 
c  determine  vertical  stress  due  to  gravity  at  zone  mid-depth 

ttempy  =  -rdn*grvyX(ypnnbr-ypnntr )X0 . 5 
ttempx  =  -rdnXgrvxX(xpnnbl-xpnnbr)XO . 5 
c  determine  complete  gravity  induced  state  of  stress 

if ( grvy . gt.0.0)call  gloadC -grp, -grp, -grp, -( trefy+ttempy+atmo ) , 

*  tyyn, txxn, tzzn, prhn, comi ) 

if ( grvx . gt . 0.0) call  gloadC -grp, -grp, -grp, -( tref x+ttempx+atmo ) , 
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225  stateCii )=exg( jj  ) 


initialize  stresses 

23 0  ep=pr ho 
p  =prho 
s(  1 ) =sxxo 
s(2)=syyo 
s( 5) =szzo 
sigxy=sxyo 
sigyz=syzo 
sigxz=sxzo 

240  call  model 

250  prhn=ep 
sxxn=s( 1 ) 
syyn=s(2) 
szzn=sC  3) 
sxyn=sigxy 
syzn=sigyz 
sxzn=sigxz 

track  invariant  shear  stress  in  eps  (zonary) 

epsn=sqrt(amaxl (atmo/1 00 . , (-2./3.X(s(l)Xs(2)+s(2)Xs(3)+s(3)Xs(l)- 
x  sigxyXX2/2.+sigyzXX2/2.+sigxzX*2/2. )))) 

calculate  total  stress  at  time  n+1 
call  zonstr 

calculate  change  in  distortional  energy  density 

call  zonzde 

prhh=(prhn+prho)/two 

calculate  change  in  internal  energy  density 
zien=ziec-£prhh+avsh)xdlrlvh+dlzdeh 

update  material  model  state  variables 
300  do  305  ii=l , nsvar(ieos) 
j j=ii*2 

305  exg( j j )=state(ii ) 

calculate  sound  speed  squared  at  time  n+1 
40C  i f ( ieos . eg . 9 ) go  to  490 

sssn=( bmodn+4 ./3 . Xgmodn)/rhoref 
go  to  500 

490  sssn=  exa*exp( -e j rlxexv ) X ( e j rl*exv*exv*( 1-exw/ ( ej rlxexv ) ) 

X  -exw/eirl) 

X  +  exbxexpC -ej r2*exv ) x( ejr2*exvXexvxt 1 -exw/ ( e j r2*exv ) J 

X  -exw/ejr2) 

X  +  exwX( exe+prhnXexv ) 

save  current  bulk  and  shear  modulii  in  zonary 

500  yldncbmodn 
shrn=gmodn 
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